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Abstract 

Within the framework of quenched lattice QCD and using 0{a) im- 
proved Wilson fermions and non-perturbative renormalisation, a high stat- 
istics computation of low moments of the unpolarised nucleon structure 
functions is given. Particular attention is paid to the chiral and continuum 
extrapolations. 
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1 Introduction 



The results of a lattice simulation of Quantum Chromodynamics (QCD) give 
in principle a direct probe of certain low energy aspects of the theory, such as 
hadronic masses and matrix elements. This is at present the only way of getting 
these quantities from QCD, without additional model-dependent assumptions. A 
useful theoretical tool in conjunction with QCD and deep inelastic scattering (or 
DIS) experiments is the operator product expansion, OPE. At leading twist the 
OPE relates moments of an experimentally measured structure function, generi- 
cally denoted by F, to certain matrix elements v n where 

I dxx n - 2 F NS (x,Q 2 ) = fE s F ., Vn [-Q^>9 S (M)j v s n (g s (M)) . (1) 

F is a function of two variables - Q 2 , the space-like momentum transfer to the 
nucleon and x, the Bjorken variable (/ is a normalisation factor). v n are the 
nucleon matrix elements of certain operators and E are the associated Wilson 
coefficients. These are perturbatively known at high energies where the coupling 
constant g becomes small and are found in a specified scheme S at scale M. 
Usually, of course, we take S = MS at scale M = \i ~ few GeV. We also 
assume that Q 2 is large enough, so that higher twist terms ie 0(1/Q 2 ) terms are 
negligible. 

As will be discussed later (section 14.3)1 . lattice computations are presently 
restricted to determining non-singlet, NS, nucleon structure functions 

F NS (x,Q 2 ) = F^x,Q 2 )-F n (x,Q 2 ), (2) 

ie the difference between proton, p, and neutron, n, results. Note in particular 
that this means that nucleon matrix elements of gluonic operators have cancelled. 

In this article we shall only be concerned with unpolarised structure functions. 
The same matrix elements v n contribute to the scattering of charged leptons 
and of neutrinos, but the weights / are different in the two cases. Thus for 
charged lepton-nucleon DIS, IN — > IX, which is mediated by a photon, we have 
F ns _^ f t,ns with n = 2,4, . . . and / 7 = 2(e 2 u - e 2 d )/2 = 1/3. For neutrino- 
nucleon charged weak current interactions for example vN — > IX, (l + N — > VX) 
or IN — > vX, (VN — > l + X) which are mediated by W + , W~ respectively, then 
we have F NS -> F^ ±]NS and f w+ = 2(-l), f w ~ = 2(+l) (neglecting the CKM 
mixing matrix) with n — 3, 5, . . .. (Alternatively setting F NS — > 2xF^ s in all cases 
one has the same matrix elements and fs as for in eq. (0) , but different Wilson 
coefficients. The additional F 3 structure functions, occuring because of parity 
non-conservation also obey eq. (0) with F NS xF 3 ,NS and again f w = =1=2 
with n = 2,4, . . ..) Similar expressions hold for the neutral currents, but with 
more complicated expressions for the fs. 
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In all cases the relevant matrix elements are given by first defining the se- 
quence of quark bilinear forms 

Of 11 - = i n ~ l q^ D^-.-D^q, q = u J d J (3) 

where D= \{D — D)- Symmetrising the indices and removing traces, gives the 
Lorentz decomposition of the proton (ie nucleon, N) matrix element of 1 

(N(p)\ \p\^-^ - Tr] \N(p)) s = 2iM )S [p^ ---p^- Tr] . (4) 
For example we have for n — 2, 



(N(p)\iq 



q\N(p)Y 



2v^ )S 



P P l" i N'l 



(5) 



and more complicated expressions hold for higher moments. Finally, the non- 
singlet, NS, matrix element is defined as 

< - V ( n )S ~ vL d)S ■ (6) 

In this article, we shall compute and in the quenched approximation 

(rif = 0), by finding the appropriate matrix elements in eq. (j3J). As will be 
seen most effort will be spent on t>2, as this is technically less complicated than 
the higher moments, and also numerically the lattice results are more precise. 
The lattice approach discretises Euclidean space-time, with lattice spacing a, in 
the path integral and simulates the resulting high-dimensional integral for the 
partition function using Monte Carlo techniques. Matrix elements can then be 
obtained from suitable ratios of correlation functions, PQI2]- Note that the lattice 
programme is rather like an experiment: careful account must be taken of error 
estimations and extrapolations. There are three limits to consider: 

• The spatial box size L$ must be large enough so that finite size effects are 
small. Currently sizes of > 2 fm seem large enough (the nucleon diameter 
is about 1.5 fm) 2 . This situation is probably more favourable for quenched 
QCD (where one drops the fermion determinant in the path integral, see 
section HJ), as due to the suppression of the pion cloud, we would expect the 
radius of the nucleon to be somewhat smaller. This indeed seems to be the 
case, see eg 0. 

• The continuum limit, a — > 0. We use 0(a) improved Wilson fermions 
(where the discretisation effects of the action and matrix elements have 
been arranged to be 0(a 2 )). For unimproved Wilson fermions, or where 
one has not succeeded in entirely 0(a) improving the matrix element we 
should extrapolate in a rather than a 2 . 



1 The nucleon states are normalised with the convention (N(p)\N(p')) — (27r) 3 2_Ejj5(p — p 1 ). 
2 Additionally all our current lattices have m ps Ls > 4 where m ps is the pseudoscalar mass. 
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• The chiral limit, when the quark mass approaches zero. There has recently 
been much activity on deriving formulae for this limit, j3J El IH1 E] 
However while most of these results are valid around the physical pion 
mass on the lattice, it is difficult to calculate quark propagators at quark 
masses much below the strange quark mass. Thus the use of these formulae 
is not straightforward. 

In addition the lattice matrix element must also be renormalised. 

Previous lattice studies gave discrepancies to the phenomenological results, 
especially for v?- In this work we want to try to narrow down the sources for this 
difference. In particular we shall present here non-perturbative, NP, results for 
the renormalisation constants (and as many previous studies used results based 
on perturbation theory compare with these other results). We also consider 0(a) 
improvement and operator mixing to enable a reliable continuum extrapolation 
to be performed. 

Compared to our previous work [2] we have improved our techniques in several 
respects: 

• We employ non-perturbatively improved Wilson fermions instead of unim- 
proved Wilson fermions. This should reduce cut-off effects. 

• Modified operators are used for v 3, t>4, which improves the numerical signal. 

• In [2 j we had simulations for a single lattice spacing only. Here we shall 
present results for three different values of the lattice spacing so that we 
can monitor lattice artefacts. 

• As mentioned before, the 1-loop perturbative renormalisation factors of [2] 
have been replaced by non-perturbatively calculated renormalisation con- 
stants. In addition we shall pay close attention to possible mixing problems 
of the operators involved. 

• We have increased the number of quark masses at each value of the lattice 
spacing in order to improve the chiral extrapolation. 

The organisation of this paper is as follows. In section El various continuum 
results for the f3 and 7 functions and for the Wilson coefficients in the MS scheme 
are collated and renormalisation group invariant quantities are introduced, while 
in section El some NS phenomenological results are discussed to compare later 
with the lattice results. Section |U describes our lattice techniques: choice of 
operators, operator mixing problems, 0(a) improvement and gives the bare, ie 
unrenormalised results. Section El discusses and compares various renormalisation 
results: one-loop perturbation theory and a tadpole improvement, together with 
the RI' — MOM non-perturbative method. In section El we discuss the results and 
give continuum and chiral extrapolations. Finally in section [7| we present our 
conclusions. 
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2 Continuum QCD results 



In this section we shall consider the RHS of eq. (jlj . Much of the functional form 
is already known: the lattice input is reduced to the computation of a single 
number (for each moment). 

The running of the coupling constant as the scale M is changed is controlled 
by the f3 function. This is defined as 



/5 5 (/(M) 



dg s {M) 



d\ogM 



bare 



-bo( 9 S ) 3 -h(g*T -W) 7 -bf(g*y - . 



,S\5 



S\7 



,S\9 



(7) 



with the bare parameters being held constant. This function is given pertur- 
batively as a power series expansion in the coupling constant. The first two 
coefficients in the expansion are universal (ie scheme independent). In the MS 
scheme where (S,M) = (MS,fi), the expansion is now known to four loops 
[TUJ [TT] . The three-loop result for quenched QCD is 



bo 



11 



bi 



102 



VMS 
u 2 



2857 



(4tt) 2 ' " a (4tt) 4 ' " a (4tt 
We may immediately integrate eq. (J7J) to obtain the solution 

r g s (M) 



(8) 



M = A exp 



1 



2b gS(My 



bog* (My 



o exp 



1 



P s (0 bo^ 



bti 



(9) 

where A s is an integration constant. 

While eq. (0) is the conventional definition of the moment of a structure 
function, for us it is more convenient to re-write it using renormalisation group 
invariant (or RGI) functions. The bare operator (or matrix element) must first 
be renormalised 



O s (M) 



J 0:bare 



(M)O, 



bare 



(10) 



giving 7, the anomalous dimension of the operator, 

d log Z a s (M) 



lo (g s (M) 



<91ogM 



bare 



doAg s ) 2 + d%Ag*Y + 4 ;2 (^) e + 



S\4 



(11) 



(The first coefficient is again scheme independent.) One may also change the 
scale and/or scheme (S, M) — » (S', M') for the operator by 



O s '(M') = Z .s'{M',M)0*{M) 



(12) 
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This leads to two anomalous dimension equations, obtained by either differenti- 
ating with respect to M' or M. Integrating these equations gives 

Z .f(M',M)= Z °f (M '' S A j (M \, (13) 



[AZg(Af)]- 1 = [2b g s (M) 2 \ 2b ° exp<M d£ 



7o(0 + d O;0 



(14) 



where we have defined 

f9 S (M) 

,/.< 

From eqs. f)12j) and (|13|). we see that we can define a RGI operator by 

O mi = AZ${M)O s {M) = Z**O bare . (15) 

Then obviously O mi is independent of the scale and scheme. The 7 function thus 
controls how the matrix element changes as the scale M is varied. Note also that 
the normalisation of O mi depends on the convention chosen for AZq(M), here 
given in eq. (fT4]) . 

As the LHS of eq. (P) is a physical quantity, the RGI form for the Wilson 
coefficient is given by 

E^SQ 2 ) = AZ s Vn {M)^E% n (^,/(M)) . (16) 

It is convenient to choose M 2 = Q 2 , as then 

E^M) = ^ s Vn {Q)-'E s F . jVn {^9 S {Q)) , (17) 

and Ep . v (1, g s (Q)) has no large numbers in it, so that a perturbative power 
series in g becomes tenable. In two schemes S and S' from eq. (|T7|) we have 

E F . Vn {l,g s {Q)) _AZi(Q) 



E% Vn (l,9 s '(Q)) A^'(Q) 



as Q — ► 00 



because in this limit g s ' = g s + ...—> 0. Hence Ep, n (l, 0) is independent of the 
scheme. With our convention for / this is 1, so that 



Ef. Vn (l,g s ) = 1 + 4 ;n;0 (/) 2 + 4 ;n;1 (/) 4 + . . . . (19) 

Practically we shall here only consider the n = 2, 3 and 4 moments. For these 
moments we have, for quenched QCD (ie rif = 0) O EH EH EH] , 

- ^{i; f; 1} -P; 3; 4}, (20) 
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Figure 1: The one, two and three loop results for [AZ^ s (n)] 1 and [AZ* 4 4 S (//)] 1 for 
quenched QCD versus fi/A 



for n = 2, 3 and 4 respectively (C(3) = 1.20206 . . .). The operator has anomalous 
dimensions given by, [TBI 113 UE] > 



dvr, 



€1 



IMS 
U v n ;2 



n = {2;3;4} 



1 J64 100 628 
(4vr) 2 IT' "9"' 45" 
1 f 23488 34450 52419141 

30375 J ' 

64486199 2200 



(4tt) 4 I 243 ' 243 
1 f 11028416 2560 



(21) 



(4vr) f 



6561 



+ 



C(3); 



+ 



26244 ' 81 
245787905651 11512 



82012500 



405 



C(3) 
C(3) 



again for n = 2, 3 and 4, respectively. Solving first eq. Q for g MS and then using 
eq. (|IH) gives the results for [AZ^ s (/i)] _1 shown in Fig. [T] Note that by loop 
expansion, we mean using the j3 and 7 function result to the appropriate order; 
we do not expand eqs. (JHJ), (|14|) any further, but solve them numerically. 

To determine a result in GeV, we shall use the ro scale here, EHI- From l20|l2T] 



we take A MS r = 0.602(48) for quenched QCD and together with the scale choice 
r = 0.5 fm = 1/(394. 6MeV) this gives for an energy of Q = /x = 2 GeV for 
example, fi/A MS ~ 8.4. The Wilson coefficient, Ep 2 s Vn can also be found and is 
shown in Fig. |21 for n = 2, 4. To obtain E^' Vn from eq. ()17j) we must simply 
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1.00 
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Figure 2: The one and two loop results for Ep£ Vn i n = 2, 4 for quenched QCD versus 
Q/A 17 *. The two loop results at 2 GeV are 1.011(1), 1.130(3) for n = 2, 4 respectively 
where the error is a reflection of the error in A MS ro. 

multiply the results from Fig. ^ with those of Fig. |21 From this latter figure we 
see that the change from the tree level result for the n = 2 moment in the Wilson 
coefficient is at most ~ 1% for Q ~ 2 GeV and is practically negligible. This is 
not so for the higher moments, when the Wilson coefficient deviates significantly 
from one. Useful values for AZ*^ S (relevant for the forthcoming lattice results) 
are given in Table |H1 in Appendix [X] 



3 Phenomenology and experimental data 

Ideally we would like to make a direct comparison between the theoretical and 
experimental result, by re- writing eq. (JTJ as, 

jf 1 dxx n - 2 F NS (x, Q 2 ) = fE™ n {Q>T • (22) 

The RHS of this equation has a clean separation between a number v^ 31 , which 
can only be obtained using a non-perturbative method (eg the lattice approach) 
and a function, Ep:' Vn (Q 2 ), which describes all the momentum behaviour of the 
moment. 

More conventional (and practical) however is to use parton densities. Usually 
phenomenological fits using parton densities are obtained from global fits (such 



S 



as MRST, [22] and CTEQ, [231 ) to the data. In this section we shall compare 
whether taking moments of the structure function gives the same answer as taking 
moments of the parton density. This could also help in estimating the error in 
the phenomenological fit. Parton densities q, q are implicitly defined by 

f 1 dxx n ~ l [q s (x, Q) + (-l)V Q)} = v^ s (Q) . (23) 
Jo 1 J 

We may relate the structure function to the parton density via a convolution. 
Defining similar but separate Mellin transformations for even and odd n by 

rl rc+ioa (Jj^l 

h n = dxx h(x) n = even/odd , h(x) = / — : x~ n h n , (24) 

(where in the inverse transformation, n in h n is analytically continued from 
even/odd integer n values to complex numbers and c is chosen so that all singu- 
larities lie to the left of the line n = c) then gives, 



F NS (x,Q 2 ) = fx C —E s (x/y, Q) 
Jo V 



o y 

x \u s (y, Q) - d s (y, Q)) + rj(lf(y, Q) - t(y, Q)) , (25) 



where r\ = +1 for even n (ie for F2' NS ) and rj = — 1 for odd n (ie for ' ns )- 
To lowest order in the coupling constant we get from eq. (fT9*j) . E s (z, Q) = 5(1 — 
z) + 0(g s ) 2 and so 3 

F NS (x, Q 2 ) = fx [(u s (x, Q) - d s (x, Q)) + ri(v?(x, Q) - 3 s \x, Q))] + ... . (26) 

The parton densities are usually determined from global fits to the data, with an 
assumed functional form, typically for MRST results like 

xq~ s (x, Q ) = A q x x " (1 - xf« (l + e q ^ + lq x) , (27) 



with parameters A q , X q , rj q , e q and 7 9 at some given reference scale Qq. For the 
MRST results given here we use the fit MRST200E, [2S1 together with the error 
analysis of >22^ at a scale of Ql = 4GeV 2 . As a comparison we also consider 
the CTEQ fit CTEQ61M, with errors calculated from [22] • (Practically, in both 
cases, we use the parton distribution calculator j2Hl to compute the moments in 
eq. (J23D.) 

Let us now briefly consider some lepton-nucleon DIS experimental results. 
While F2 ,p (x, Q 2 ) is well known, experiments with deuterium to find ^'"(x, Q 2 ) 
are much more difficult, due to target nuclear effects. We shall use here the 
results from [27] which employ both proton and neutron (deuterium) targets in 
the same experiment, which thus minimizes systematic errors. In |2E] this has 
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Figure 3: World experimental data for FJ' (x, Qq), 28 , at = 4GeV^ in the form 
of bins, plotted against x using a linear scale. Errors in the bins are also shown. The 
dot-dashed line is a rough estimate, obtained by a linear extrapolation of the last bin 
to nought (at x = 1). 



been combined with the world data, |2H|, and is shown in Fig. El in the form of a 
series of bins at different x values. (Naively, if there were no QCD interactions, 
the parton model would give a delta-function distribution at x = 1/3. This 
distribution has been considerably washed out here though.) There is a paucity 
of data for larger x. However _F 2 y ' iVS is dropping rapidly to zero, so any error here 
will not affect the low moments. As shown in the figure, we have simply made 
a linear extrapolation to x — 1, to estimate this region. To find the moments 
for eq. (fl)l. we simply need to find the area under the bins weighted with the 
appropriate x-moment, ie 



J ® J bins.b 

j 0.164(4)(1) n = 2 
~ i 0.0289(10)(10) n = A ' 



(2c 



where the first error is from the data and the second error is the effect of dropping 
the estimated last bin. (As expected, we see that higher moments are more 
sensitive to this bin.) 



3 There are many references to the relationship between structure functions and parton 
densities. See for example |24|. 
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'World' 


MRST 


CTEQ 


< s (Qo) 

„H3T 


0.161(4) 
0.226(14) 


0.157(9) 
0.220(18) 


0.155(17) 
0.217(27) 


< s (Qo) 




0.0565(26) 
0.0972(95) 


0.0551(51) 
0.095(12) 


< s (Qo) 

U4 


0.0241(13) 
0.0480(58) 


0.0231(11) 
0.0460(54) 


0.0230(23) 
0.0458(67) 



Table 1: Values of v„ s (Q ) at Q = 2 GeV. The Wilson coefficients (needed for 
World') have been computed from eq. Q19|). using A MS | ra/= 4 = 250(50)MeV, giving 
Ef 2 -v2 = 1-018(3), Ef 2 -,v4 = 1-200(30), where the error is a reflection of the error in 
A MS r . Similarly, to convert to the RGI form, [AZ^f(2 GeV)] _1 | n/=4 from eq. (fT5|> 
again uses Kf=± giving 0.713(40), 0.581(50) and 0.502(54) for n = 2, 3 and 4 respec- 
tively. 

In Table [U we give estimates of v^ IS (Qo), v^ s (Qo), using estimates of the 
Wilson coefficients and [AZff s (2 GeV)] -1 given in the figure caption. (These 
numbers are similar to the quenched results, as can be seen from the caption 
of Fig. Eland Table 0) We find that there is good agreement between the two 
methods, with the lowest moment from MRST or CTEQ being slightly smaller 
than the experimental result. Thus these global fits describe the (low) moment 
data well 4 . In future for definiteness we use the MRST results. 



4 The Lattice Approach 

Euclidean lattice operators 5 are defined by 



O 



q ■ 



(29) 



where q is taken to be either a u or d quark and V is an arbitrary product of 
Dirac gamma matrices. (The q index will be usually suppressed.) We have used 
the lattice definitions 



D„ q(x) 



1 r 

2a 
1 



UJx)q(x + a/2) — U*Jx — ajX)q(x — a/2) 



q(x) Dn = — q(x + afyU^x) — q(x — afyU^x — a/2) 



(30) 



4 Note that a recent analysis, jSH), gives similar results of v^ s (Q ) w 0.159 and v^f s (Qo) 
0.0245. 

5 Our Euclidean conventions are described in 
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and D= \(D — D). For the operators corresponding to eq. Q, we shall only 
need V — 7. However for the discussion on mixing under renormalisation, we 
shall also use Y = 775 and o"7 5 . 

4.1 Choice of lattice operators 

We now take the simplifying choices of two momenta p = or p = (pi, 0, 0) = pi 
with p\ being the lowest non-zero momentum possible on a periodic lattice ie 
Pi = 2tt/Ls where the number of sites in each spatial direction is Ls/a. We take 
our lattice operators as 6 , 

CP 

^{14} > 

CP - I ( CP + CP 4- CP \ 

^{44} 3 V {11} + ^{22} + ^{33} J ' 
C{441} _ 5 (^{221} + ^{331}) ' 

^{1144} + ^{2233} ~~ \ (^{1133} + ^{1122} + ^{2244} + ^{3344}) ' (^1) 

Of course, there are other possibilities. However these will all require non-zero 
momenta in two directions or suffer from more severe mixing problems. As we 
shall see, even a non-zero momentum in one direction leads to a strong degrada- 
tion of the signal and with two non-zero momenta very little signal is observed, 

The transformation properties under the hypercubic group H(4) are given 
in Table |2l Note in particular that the 'off-diagonal' (0 V2a ) and 'diagonal' 
operators (0 V2b ) for t>2 belong to different representations, in distinction to the 
continuum operator. Thus we expect that although these operators should have 
different lattice artefacts and renormalisation factors, in the continuum limit both 
should lead to the same result - potentially a useful check. 

4.2 Mixing of lattice operators 

A given operator of engineering dimension do can mix with operators (with the 
same quantum numbers) of lower dimension, the same dimension or higher di- 
mension. 0(a) improvement involves mixing with one dimension higher operators 
(irrelevant operators) where the choice of improvement coefficients is convention- 
ally treated separately to mixing with operators of dimension < do (relevant 
operators). We shall follow this practice here. 

e O V2a , V2b are the same operators as we used previously in 0, while there are some mod- 
ifications to V3 and V4 . For V3 we have effectively made the proper rotation (ie preserving 
75) of 4 <-> 1 (and 2 ^ 3) of the operator in 0. This means that the measured ratio in eq. (|39(l 
is —ipiE^v^ rather than — p\vz and hence the signal is better by a factor 0(Ep 1 /p\). For 
V4 we have symmetrised over the 2, 3 indices of the operator given in |2], which makes the 
measurement of the ratio for the new operator numerically a little less noisy. 
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l'2a 



o 



''■2b 



o 



"3 



a 



12 



Operator 


(rl",C) 


o V2a 


(rf,+) 
(rf»,+) 


o 


(n (2, ,+) 



Table 2: Transformation of the various operators under the hypercubic group H(4), 
|33j . where I represents the dimension of the representation rj. , k labels different 
representations of the same dimension and C is the charge conjugation parity of the 
operator. 



4.2.1 Operator mixing with additional relevant operators 

While for V2 there is no mixing with relevant operators, unfortunately for O v:j 
and Vi there are other relevant operators transforming identically under H{4) 
which can thus mix with the original operator, [33]. Specifically we have 7 

• Operators mixing with V3 : 

^3 ~ ^((144)) 2 ^((133)) + U ((122))) ' 

^v 3 ~ ^||432|| OLy |432| — ZLy ||342|| > l 0Z J 

where using the notation of [33] we have denned: 



20 T + 20 T 

(33) 



\v\v2v3\ 


= o T 

V\V 2 V3 


^1/11/31/2 




+ o r 

1 ^1/31/21/1 


\Wiu 2 iy 3 \\ 


= o v 




+ o T 






= o r 


+ o r 

1 "1 




1/3 "2^1 



7 Using the operators given in [5| we would have 

0">i — f)T - I I f) 7 4- fl 7 ^ O" 12 — JO 115 

^v 3 —^((411)) 2 ^<(422>) + ^((433))^ ' ^v s — Z ^||213|| 
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Operators mixing with O 



(nmi 

V 4 



CV — cp — rp — rp _i_ on 7 -i- 9C9 7 

^1144 ^4114 ^1441 ^4411 ' ZLy 1414 T" -^4141 

^3322 + 202323 + 203232 
+ \ ( + <?7l33 + <%13 + <%31 + O? 



"^2233 ~ ^3223 — ^2332 ~ ^3322 + 202323 + 2£ 

~2 ( + ^1133 + ^3113 + ^1331 + ^3311 ~~ 201313 

+ ^1122 + ^2112 + ^1221 + ^2211 ~~ 20 7 2 12 

+ 02244 "I" ^4224 "I" ^2442 "I" ^4422 — 202424 
/y\*y //^\ r y c\ /y~\ 'Y 



ZLy 3131 

0/-O7 
ZLy 2121 



- 2(9 7 3 4 3 J 



«4 



0^3 

^4 



I, 

2'' 



o/n7 
ZLy 4242 

— 344a 1 ~ 4433 ^^3434 
^1432 "T ^3412 "I" ^2143 — ^4123 — ^2341 "I" ^4321 

fills , /rpns i /nTTB /o^ 5 _l ro^T 5 _l roTT 5 

^1324 1* ^2314 1* ^1423 ^2413 ^3142 ~>~ ^4132 "1 ^324 

1 A1775 /n77S /'Ol">' 5 _l_ ^OTT 5 " " """"" 

"t" Ly 1243 ^4213 ^1342 ' ^431: 



-)7 , /-o7 
7 2442 1 ^4422 

+ 03344 + ^4334 + ^3443 + ^4433 

l"r L/ 1234 — ^3214 _ 



1 + ^2413 
^^^2314 



/O "^ 5 _l_ ^0°"75 
^2431 ^1324 

_ 0/90-75 I 9/-9075 

ZL/ 2341 ^ ZLy 1423 



^3124 

/O "^ 5 _l_ <*T T7s 
^3412 1 ^3421 

ZL/ 1432J • 



)775 1 ^)775 
4312 1 ^2134 

^1342 



^4231 ) 

/■nils 1 /--m 5 

^2431 1" ^3421 

/'O "'*' 5 _l_ /O ct 75 
" ^1234 1 ^1243 



(34) 



O™ 3 is an operator with one lower dimension (and different chiral prop- 
erties) than the other 0„ 4 operators. It is also possible that four-fermion 
operators can mix: we shall not consider this here though. 

In Appendix |B] we illustrate, by an example for 0„ 4 , how 0™ 1 and 0™ 2 can 
be derived. The other mixing operators follow by similar considerations. 

While this list contains all operators allowed by group theoretical arguments, 
not all the operators contribute, see sectional 

4.2.2 Operator Improvement 

As we are using 0(a) improved fermions, to achieve the elimination of 0(a) terms 
in matrix elements, the corresponding operators must also include additional 
irrelevant terms, with coefficients appropriately chosen. Presently only for the 
lowest moment (ie v 2 ) are these extra operators known, [33]. In this case we 
should replace the operators 7 in v 2a or v 2 b by 8 



7 

jJtV 



i=l 



1 + arriqCo) qjf, D v q 



+ac\ 



+ac 3 



i J2 ( l (J ^ D[vDx\ q 
x 



L A 



qcr^x D v q 



+ ac 2 



-q D{^Du} q 



(35) 



8 0ften, alternatively, the improved operator is re- written as (l + am q bo)(0 + ^2 i=1 ac^O^). 
In the remainder of this article it should hopefully be clear from context whether we are referring 
to the operator or its improved version. 
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where m q is the bare lattice quark mass, related to the hopping parameter k by 

am * = \ {I - £) • (36) 

So we see that there are potentially four additional improvement operators, de- 
fined in eq. ()35|) together with four unknown improvement coefficients, Ci(go), for 
i = 0, . . . , 3, which are functions of the gauge coupling constant g . The i = 3 
operator only contributes for non-forward matrix elements, which are not consid- 
ered here and so this term may be dropped. Also ultimately Co will not concern 
us as we are only interested in the result in the chiral limit. For on-shell matrix 
elements the equation of motion may be used to eliminate one of the improvement 
terms, for example we can choose Cq and C\ as linear functions of c 2 . First order 
perturbation theory gives the relation between the improvement coefficients for 
v 2b of 



c = l-c 2 + ^4(17.20377 - 8.69045c 2 ) + O^, 

ci = c 2 + O(g 2 ), (37) 

{Cp = 4/3) and a similar expression for v 2a , |H11- With these values of the 
improvement coefficients, 0(a) corrections to the v 2 b matrix element have been 
eliminated (at least to lowest order perturbatively). We see that we cannot 
determine the O(gg) term of C\ because it is the coefficient of an operator that 
vanishes at tree level. Non-perturbatively the improvement coefficients have not 
yet been determined. However one might suspect, that choosing c 2 = also gives 
in this case a small c\ coefficient. 

For higher moments, the bases for the improved operators become increasingly 
cumbersome, as not only would we expect more possible irrelevant operators built 
from the original operator together with an additional covariant derivative, but 
also four-fermion operators may play a role. This does not necessarily detract 
from the original operator though, because we can always attempt to make a con- 
tinuum extrapolation in a rather than a 2 if we cannot motivate why the irrelevant 
matrix elements are small. 

4.3 Determining the matrix element 

Matrix elements are determined from the ratio R of three point to two point 
correlation functions, 

R(t, r; Pi O) = y 74 ) 77 a , (38) 

where Ci. x At\ p) is the unpolarised nucleon two point function with a source 
at time and sink at time t, while the also unpolarised three point function 
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Figure 4: The 3-point quark correlation function for a baryon. 
Ci ^ ^(t,T;p;0) has an operator insertion at time r. If we consider a region 
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<C t <C t < \L T (ie well inside the nucleon branch of the propagator) where 
Lt is the temporal extension of the lattice, then the transfer matrix formalism 
upon projecting out the ground state nucleon leads to 

R(t,r;p;0 V2a ) = ipiv 2a , 

p2 I 1^2 



R(t } r;p;O v J = P —^-v 2b 

bp 



R(t,T;pi,0 V3 ) = -ipiE^v 3 , 

R(t,r;px;0 V4 ) = p\E^v±, (39) 

for the bare matrix elements v n . For some more details see eg jSHI2j- in general 
we would expect that the best signals with the smallest noise are seen for zero 
momentum. We take (and have numerically checked) that the standard dispersion 
relation E~ = m 2 N + p 2 holds for p '= p\. 

The nucleon three-point correlation function is depicted in Fig. 0] While the 
two-point correlation function consists of one diagram, for the three-point corre- 
lation function we have two diagrams - a 'quark line connected' contribution, and 
a 'quark line disconnected' contribution, left and right diagrams in Fig. ^respec- 
tively. This is not the usual field theoretic splitting of the Green's function into 
connected and disconnected diagrams. As quarks can travel backwards in time as 
well as forwards, we would expect that the quark line connected term would also 
give a contribution to the anti-quark parton density defined in eq. (j2Sl)- As quark 
line disconnected diagrams can, by definition, only interact with the hadron via 
the exchange of gluons then the numerical results suffer from large short distance 
(ie ultra-violet) fluctuations. So a very large number of configurations is required, 
which is very expensive in computer time. We have not computed this term here. 
To cancel any effects of these disconnected terms, if flavour SUp(2) symmetry 
is a good symmetry, it is sufficient to consider NS matrix elements, ie the u 
quark matrix element minus the d quark matrix element. In Appendix [CJ for 
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completeness, we give explicit expressions for the relevant three-point correlation 
functions. 

A further class of disconnected terms are given by gluon matrix elements. 
Again on the lattice these are difficult to compute, see [35] , but again they cancel 
upon considering non-singlet matrix elements. 

All these gluonic or sea-quark effects are concentrated at small x, and thus 
for higher moments, ie n = 3 or 4, are naturally suppressed anyway. Thus 
disconnected contributions may be less significant, so that the computation of 
singlet matrix elements is then more reliable and we can consider just a u or d 
operator matrix element. Although the quenched approximation does not handle 
the sea-quarks correctly, we might also expect for these higher moments that 
quenching has less effect. But these statements are hard to quantify, and as this 
is all less likely to be the case for the lowest moment anyway, we shall consider 
mainly the non-singlet results here. 

4.4 Raw results for lattice matrix elements 

We now discuss our raw numerical results for the lattice operators and the nu- 
merical significance of the additional improvement operators and/or additional 
relevant operators to the nucleon matrix element. Since our original publica- 
tion [2] which employed unimproved Wilson fermions at (3 = 6.0, we have used 
quenched configurations with 0(a) improved fermions: thus for the action we 
take the standard Wilson gluon action, while for the fermion propagator the 
standard Wilson fermion is used together with a 'clover' improvement term. The 
(non-perturbative) coefficient c sw of the improvement term was taken from [36J. 
This means that on-shell quantities, such as masses, have only 0(a 2 ) discretisa- 
tion effects. Simulations were performed at three (3 = Q/g^ values, f3 = 6.0, 6.2, 
6.4. At each f3 value, four or more quark masses (degenerate in u and d) were 
used, at each mass a statistic of several hundred configurations was generated. 
Antiperiodic fermion boundary conditions were taken in the time direction and 
periodic in the remaining spatial directions. Further details of our runs are given 
in Appendix [01 in Table [7| 

To improve the overlap with the nucleon we employed Jacobi smearing, and 
used a non-relativistic (NR) projection of the nucleon. Jacobi smearing is de- 
scribed for example in |31j . where a hopping parameter (k s ) expansion (of order 
n s ) of a Wilson fermion operator restricted to a time plane smears out the orig- 
inal point quark source. For (3 = 6.0, 6.2 and 6.4 we use (n s ,n s ) = (0.21,50), 
(0.21, 100) and (0.21, 150) respectively giving a root-mean-square radius of about 
0.4 fm, a reasonable fraction of the nucleon radius ~ 0.8 fm. The NR projection, 
where in our Dirac gamma matrix representation only the upper two components 
of the spinors are used, is briefly described in [37] and Appendix IU1 

The nucleon sink positions were chosen as t — 13a, 17a and 23a for f3 = 6.0, 
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Figure 5: V2b versus r/a from the ratios R( 17a, r; 0; Oy 2b ), R(17a,T;pi;0 V2b ), eq. l|39j) 
left and middle pictures respectively and via from i?(17a, r;p*i; V2a ), right picture 
for (empty circles) and NS (ie 0( u ) - 0W) (filled circles) for (3 = 6.2 at 

k = 0.1344. The chosen fit intervals are denoted by vertical dotted lines. 



6.2 and 6.4 respectively and the fit ranges in r were taken as [4a, 9a], [6a, 11a] 
and [7a, 16a]. These fit range values were not so critical, but allowed a splitting of 
the range [0, t] into roughly three equal parts each piece being roughly the same 
physical size. 

In Appendix |Dl the data is presented in tables giving separately the u and d 
contributions. As discussed previously though, as the disconnected diagrams have 
not been computed, the most physically significant result is for the non-singlet 
matrix elements. Thus we have also repeated the data analysis directly for these 
matrix elements. This allows, in particular, a better estimation of the error. 
These numbers are also given in the tables in Appendix [Dl As an example of the 
typical ratios obtained, in Fig. El we show R(17a, r; 0; V2b ), R(17a,T]pi]0 V2b ), 
R(17a, r;pi, V2a ), left to right pictures respectively, for (3 = 6.2 and k = 0.1344. 
We seek a plateau in the region <C r/a <C 17. The region chosen is denoted 
by vertically dashed lines. Clearly the operator corresponding to t>2& delivers a 
better signal for the ratio than the operator for u 2a , although even for the V2b 
operator we see that it is better to choose zero momentum rather than non-zero 
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momentum. 

In this section, as mentioned before, we wish to merely estimate the numerical 
significance of extra operators, as described in section EOl To this end, as mixing 
coefficients and improvement coefficients are not so well known, we make a series 
of plots either comparing the ratio of the additional improvement operators (ie 
the v n constructed using the i n eqs. and (JHHj)) to the original operator, 

< n = ^ , (40) 

or compare directly the v™ 1 derived from eqs. ()32|) . ()34|) to the original operator, 
v n . (It will be seen in section El that the ratios r l Vn have a physical significance.) 

These are all plotted against the square of the pseudoscalar meson mass. 
As (am ps ) 2 is proportional to the quark mass for small quark mass, using this 
quantity avoided the necessity of first determining the critical hopping parameter 
k c in the quark mass definition, eq. (}36|) . We have thus used the (dimensionless) 
extrapolation parameter (r /a) 2 x (am ps ) 2 = (rom ps ) 2 , with r /a being given by 
the Pade formula in j^HJ - Using tq = 0.5 fm this enables us to get an idea of how 
close our simulation points are to the chiral limit in physical units. 



4.4.1 v 2 

In Fig.Elwe show ar\ a , ar\ a plotted against (r m ps ) 2 , together with a linear chiral 
extrapolation. Also plotted is the approximate value of the pseudoscalar mass 
corresponding to degenerate u/d or s quark masses. So we see that our simulation 
runs over the range from about two to three times the strange quark mass to a 
little under the strange quark mass. It is also noticeable that we are a long way 
from simulating with a light u/d quark mass - in fact within our errors, there 
is no difference between linearly extrapolating to the chiral limit or to the pion 
mass. Nevertheless we see that the effect of any extra improvement operators is 
likely to be small, of the order of a few percent. 

The above result is for the off-diagonal operator, which needs a non-zero 
momentum in its evaluation. This figure is to be compared with the result using 
the diagonal operator, which advantageously may use zero momentum. In Fig. [7| 
we show this result. In comparison with the previous picture, the errors are 
considerably reduced, as p = is used (see also Fig. EJ). Indeed using p = p\ (see 
Appendix \D\ for the numbers) we see that the errors grow again, although they 
are never as large as for the off-diagonal operator. Again, the effect of any extra 
improvement operators is small. 

Thus all the numerical results and linearly chirally extrapolated results for 
ar l V2 look small giving av\ < 5% of v 2 and some indeed are consistent with zero. 
From eq. (|37|). we see that if we choose c 2 = then C\ = 0(<7q), which is likely 
to be small. So these results lead numerically to a small additional improvement 
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Figure 6: The ratios ar* , ar 2 2 plotted against (rom ps ) 2 . The filled symbols at non- 
zero (rom p <;) 2 are for arl 2a , while the empty symbols denote cir 2 2a . Also shown is a 
linear chiral extrapolation (dashed line and full line for arl 2a , ar^ 2a respectively). The 
result in the chiral limit (rom ps ) 2 = is also indicated, again using full, empty symbols 
for arl 2a , ar% 2a respectively. The small horizontal lines around the chiral limit represent 
the perturbative estimate, see sectional with again a dashed line corresponding to arj 
and a full line to ar% 2a respectively. The rough values of the (hypothetical) pseudoscalar 
mass composed from u/d or s quark masses, evaluated from m w and mx respectively 
are shown as dashed vertical lines. 



term and we conclude that the effect of these terms in t>2 can be assumed to be 
negligible. 



4.4.2 t> 3 and v 4 

We now present our results for the higher moments. In Fig. |H] we show v™ 1 , 
v™ 2 for j3 = 6.0, together with a linear chiral extrapolation. Also shown, for 
comparison, is the operator v^. As expected while the magnitude of the noise has 
increased in comparison with t>2 and scatters more, an acceptable signal is still 
seen. We find a clear separation between v% and the mixing operators (indeed 
they are consistent with w 0). For higher (3 values, the data fluctuates more and 
it becomes more difficult to disentangle the results. 
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Figure 8: Mixing terms, v™ 1 , v™ 2 for (3 = 6.0 (filled circles). Also shown is v% (open 
circles) and linear chiral extrapolations. The same notation as for Fig. |HJ 
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Figure 9: t> 4 for j3 = 6.0 together with v™ 1 , i;™ 2 , v™ 3 /a. The notation is the same as 
for Fig. El 



In Fig.Elwe plot v™ 1 , v™ 2 , v™ 3 /a, together with t> 4 . While there is a reasonable 
signal for t> 4 , the improvement terms fluctuate a lot (again becoming worse the 
higher the value of (3 is). Indeed even finding a plateau for the ratio R(t, r;p; O) 
in eq. (|55|) becomes problematical. It would appear that while, numerically, v™ 1 , 
v™ 2 are small in comparison with t> 4 the situation is less clear cut for v™ 3 /a 
(although it is only including the heaviest quark mass point that leads to a non- 
zero result in the chiral limit). Nevertheless as most of the results for the mixing 
terms are much smaller than t> 4 we shall ignore any effects from them. 



5 Operator renormalisation 

A lattice operator (or matrix element) must, in general, be renormalised. Again, 
we shall discuss mixing and 0(a) operator improvement separately. 
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5.1 Operator mixing and renormalisation 

For operator mixing we can generally write 

Of (M) = £ Z 0i0 ..f Ar (M, a)O^(a) , (41) 

j 

working in a scheme S at scale M. So we must determine a matrix of renormal- 
isation constants. While for those mixing operators of the same dimension as 
the original operator, low order perturbation results can be calculated (at least 
in principle), for lower dimensional operators this is not reliable: the renormali- 
sation constant is proportional to a positive power oil /a and non-perturbative 
terms may contribute. For the higher moments, both cases occur. For V3 both 
the mi and 7712 operators have the same dimension as the original operator, and 
for t?4, mi and 777.2 have the same dimension, while is dimension one lower, see 
eqs. (j22D, (HI). 

5.1.1 Renormalisation and relevant operator mixing 

In this section we want to make a few comments on the operator mixings seen for 
the operators we use in this paper, and contrast them with the mixing problems 
for the operators which we rejected. 

In the continuum, symmetry under 0(4) in the Euclidean case (or the Lorentz 
group in the Minkowski case) imposes strong constraints on which operators 
can have unpolarised forward nucleon matrix elements. One way of stating the 
rule is to say that the only operators with non-zero unpolarised forward matrix 
elements are those which have J p = + in the nucleon's rest frame. The quantum 
numbers J p = + only occur in 0(4) representations of the form JJ^'^ 1 in the 
notation of 39j . Considering the 0(4) classification of the operators is relevant, 
because although bare lattice matrix elements only respect the symmetries of the 
hypercubic group i?(4), renormalised quantities should respect the full continuum 
symmetry group. 

If we look at the mixing operators for V3, as listed in eq. (}3*2*j) . we see that 
under the hypercubic group they transform exactly the same way as V3, but under 

0(4) they both transform according to the representation [A2'2<' ) which means 
that by the continuum symmetries their renormalised forward matrix elements 
must be zero, 

v^ s = . (42) 

Similarly the mixing operators for t> 4 , listed in eq. (JMJ), all transform according 
to the [/(°' 2 ) representation of 0(4), so 

v^' s = . (43) 
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If the renormalised mixing operators all have the value zero, we can express all 
the bare lattice matrix elements in terms of a single renormalised matrix element. 
Again taking the V3 system as our exemplar, 



Z- 1 ) mi v s 3 , (44) 

/ v. ' ;v 3 



(temporarily dropping, for clarity, other arguments of the renormalisation con- 
stants). Since all the lattice matrix elements are multiples of the physically 
interesting matrix element v% we have a choice in how we calculate the renor- 
malised matrix element. We could add up all the terms in eq. ()41)1 . assuming 
that we know the mixing Zs well enough. Or we could equally well just invert 
eq. (jUj) and calculate t>f from t> 3 alone, 

v " = Tz^~ ■ (45) 



If we calculate the renormalised matrix element from t> 3 alone, we should really 
renormalise by dividing by (Z~ x ) v ^ vz instead of multiplying by Z V3 . V3 . In practice 
the difference is minor, 



(z-i 



' -"Zw-E ^f 2 **'** - (46) 



The difference involves the product of two off-diagonal Zs, and so is 0(<7q) i n P er ~ 
turbation theory, see section lo.2.11 and is still likely to be tiny non-perturbatively. 
Note that eq. ()46|) tells us that mixing with lower dimension operators is no more 
dangerous than mixing with operators of the same dimension, in a case like this 
where the lower-dimensional operator has zero renormalised matrix element. This 
is because the product Z miZ v mi is dimensionless, even when v™ 1 is a lower 
dimensional operator. 

We conclude that for the operators we use in this paper, mixing is relatively 
benign, because continuum symmetry says that the renormalised mixing matrix 
elements are zero. 

Finally we want to contrast this with an example where this is not so, to 
show the importance of choosing the operators carefully. We could have tried to 
measure V3 with the operator 

C{ 444 } — 0(411} — C| 422 | — 0| 433 | . (47) 

This however would have much worse mixing problems because it mixes with 
operators which are allowed to have nucleon matrix elements in the continuum, 
for example 

0J 4 44 } + 0J 4U} + 0J A22} + Oj m , q^q or qD^q. (48) 
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Now we would not have the option of using eq. (JSJ), we would have to use eq. (J3TJ) 
and to get a reliable answer we would need to know the off-diagonal elements of 
Z accurately (especially those for the operators of lower dimension) . This is why 
we rejected the operator eq. (|47jl. and why we think that the mixing operators of 
subsection 14.2.11 are not a problem. 



5.1.2 Renormalisation and operator improvement 

For 0(a) improved operators, we know that the operator depends linearly on the 
improvement coefficients q as shown in eq. Our experience from extensive 

perturbative calculations in [33] leads us to expect a result of the form 



Agf = ^A^ + 0(a°) , (49) 
a 

where A 1 ^ is the lattice amputated three-point Green's function, O is the relevant 
operator and O l are the improvement or irrelevant operators in eq. (|35p. (Note 
that there is nothing forbidding the improvement terms mixing with O/a.) The 
only case that does not mix is the mass improvement term am q O, 

arriqA^ = 0(a) . (50) 

We can calculate the q dependence of the renormalisation constant by requiring 
that the renormalised A^ should be independent of the Cj at leading order in a 
and thus we have (temporarily suppressing additional M, a arguments) 

AS = Z ° ; H {Ci}) [A^ + E^Ag°- + 0(a) 

S 



Zo;LXr({ c i}) 



(l+5>Ci)A£f + 0(c 



7 S 

^#V + 0(a), (51) 



and so we can write 



O s (M) = Z s L4r (M, a, { Cl })0(a, {q}, c ) = %^|^ 0(a; {c,}, c ) , 

(52) 

where the coefficients Q = Q(go) = O(g^) have to be determined. Again in [32] 
we have given these coefficients for vi a and v-s, in lowest order perturbation theory 
using eq. (|4"5j) . (Requiring 0(a) improvement also determines the q as given for 
example in eq. (|37JI .) Numerically these coefficients turn out to be quite small. In 
addition to the perturbative Q's we can also find non-perturbative Q's by looking 
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at our measured nucleon matrix elements and requiring that the improved result 
and the unimproved agree to leading order in a, ie from eq. (|52|) we have 



1 + J2i^0 C iC 

giving, cf eq. 



v n + 0(a) , (53) 



< = -Ci«n + O(a ), (54) 
a 

or 

G = <J _+0(a), (55) 

CLTYlq — U 

(see eq. f)4U|)). We would expect Q calculated in this way to agree up to ambi- 
guities of 0(a). In Fig. Eland Fig. Qwe compare the values for Q obtained this 
way with the 1-loop perturbative results from j33] (shown in the figures as short 
horizontal lines about the chiral limit). Although, not unexpectedly, there are 
differences to the perturbative result, the ratios remain small. From eq. (|37j) we 
may choose C2 = and so c\ = O(g^) is also small. Hence the change in the 
denominator of eq. (|52|) from 1 is small and so does not change the renormalisa- 
tion constant perceptibly. Thus we shall ignore any small effects here (just as in 
section [Ol where our conclusion was that we could numerically drop the 0(a) 
improvement terms). 



5.2 Determining the renormalisation constants 

To define the renormalisation constants a renormalisation procedure must be 
prescribed. Often the renormalisation constants are defined first in a MOM 
scheme by computing the (Landau) gauge fixed two-quark Green's function with 
one operator insertion and setting 

[A£H p2 = M 2 = A£H^ , (56) 

where Aq om is the (renormalised) amputated three-point Green's function for the 
operator O. The RHS of this equation is the tree level value (or Born approxi- 
mation) of the amputated function. This definition may be used for perturbative 
computations, see eg (30] > where [. . .] in eq. (|56|) means that, in a given basis, 
we drop terms not proportional to the Born term. This may be modified, using 
a trace condition for the definition of the renormalisation constants, to give the 
alternative RI' — MOM scheme which is also suitable for non-perturbative 
calculations. A discussion and some results for this non-perturbative method will 
be given in section 15.2.31 

The resulting Z . M ™ may, if wished, be converted to another scheme, ie to the 
MS scheme using eq. ()12|) where Z .. ™f OM (n, /i p ) is (perturbatively) calculable. 



26 



The MS scheme is particularly convenient, as the renormalisation constants are 
independent of the gauge fixing condition chosen. 

Unfortunately, the definition given in eq. (JHHj) has its limitations: O™ 2 , O™ 1 
(9™ 2 and O^f all have vanishing Born matrix elements between quark states 
(as they involve commutators of covariant derivatives). So in order to be able 
to compute the renormalisation constants we would have to consider more gen- 
eral Green's functions (eg quark-gluon) . At present we must simply ignore this 
problem though. 

After these more general remarks, we shall now give various results for the 
renormalisation constants for (one-loop) perturbation theory, TRB perturbation 
theory and finally a non-perturbative determination of the relevant constants. 



5.2.1 Perturbation Theory 

One loop perturbation theory 9 yields 10 

Zo^LiM, a) = 8 0i , 0j + g 2 Q \d 0i ;O r ,o ln(aAf) - B s 0i . Q (c sw )\ + 0(g$) . (57) 



In the MS scheme [M = //) we have, 

BZ s (c„„) = — ^ r 1.27959 - 3.87297c,,,, - 0.67826c! 



v 2a \"sv)j (47r) 2 V v.^.~~.~sw "-"'^"^ > 

BT{c sw ) = (2.56184 - 3.96980c sw - 1.03973cM , (5* 



(4tt)< 

with Cp = 4/3. The calculations have been extended by S. Capitani [33] to now 
include B^ 3 S , B^ s . For V3 , off-diagonal elements in eq. (|57jl have also been 
computed 

B^(c sw ) = (-12.12740 - 2.92169c sw - 0.98166^) , 

r< 

B™,mi (c sw ) = (-0.36848 - 0.032760c sw + 0.029137c* w , . 



V3;v. 



Bffi(c sw ) = (-14.85157 - 2.15228c sto - 1.70741^ , 

B%fi. V3 {c sw ) = (-3.30605 + 0.33335c sw - 0.37050cL) • (59) 



For Bt w we have 11 



U4 



B™ s (c sw ) = j^- (-25.50303 - 2.41788c sw - 1.12826c*,,) . (60) 



9 For a general review of lattice perturbation theory, see 021 ■ 
10 For diagonal elements we write Bq.. . = B^. and do^Ovfi = ^e>;;0- 

11 The number for c sw = was incorrectly found in We shall use the result of [23| m 

the following. 
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The lowest order anomalous dimension coefficient, being universal, is given for 
the three moments by eq. (}2*T]) and in our basis the coefficients d Va . v ™i. and 
d v ™i. V3 . = vanish, while d v p. = -28/(9(47r) 2 ). 

For consistency at this order in perturbation theory, we take c sw = 1 (the tree 
level value). Setting fi = l/a to avoid large logarithmic factors gives, for example 
at j3 = 6.0, the results 



y MS 


1.0276, 




y MS _ 


1.0207, 




y MS 


1.1354, 


Z mi MS = 0.003142 , 


y MS 
V4;LAT 


1.2453. 





While we see that for v-i first order perturbation theory changes the tree level 
result (= 1) very little, there are perceptible differences for the higher moments. 
Note also that the mixing renormalisation constant for v 3 is very small in compar- 
ison to the diagonal renormalisation constant, Z In addition, although the 
mixing operator signal is rather noisy, v™ 1 v 3 as we have seen in section 14.4.21 

Thus assuming that for a non-perturbative evaluation Z V3 .^, ^> Z m-i ( clS IS 

the case for the perturbative result), we can ignore the effects of the mixing term 
in the future. 

We also use the three loop result from Table El for AZ* s (/i) to find the RGI 
factor Z£f, eq. (|T5j) . 



5.2.2 TRB perturbation theory 

To improve the perturbative renormalisation results of the last section, we shall 
apply tadpole- improved renormalisation-group-improved boosted perturbation 
theory or TRB-PT, [13], which we shall now describe. The renormalised op- 
erator is given by 

CF a {p) = Z ™{ji,a)0{a). (62) 

As in eq. ffTTj) . we may define a 7 function either in the MS scheme or what 
we shall formally call here the LAT scheme. Additionally as expansions in the 
bare coupling constant seem to be badly convergent, we choose to expand in the 
boosted coupling constant and thus we have 

= ltf(9a)=d o . g 2 n + d^ 1 g i D + ... , (63) 
where 

0& = 4. 4 = (lTrU n ), (64) 
u 



togZ ^(n,a) 



d log a 
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o 


JLAT 
a O;l 


V2a 


-152.14/(4tt) 4 


V2b 


-176.31/(4tt) 4 


V3 


91.828/(4tt) 4 


Va 


382.32/(4tt) 4 



Table 3: Values of dffil for 0{a) improved fermions from eq. 1)67(1 . 



and U is the product of links around an elementary plaquette. Expanding uq 
we have w — 



X 9oPi + O(go) where p 1 = 1/3. 



From eq. (JI3J), ie integrating eq. QJ for (5, M) = (MS,/j) and (LAT, a -1 ), 
eq. (jHSI), gives 

(65) 



MA"/ 



and thus from eq. ((62)) . the RGI quantity may be written as 
C M = AZ|P(/i)(9 w (/i) = AZ^(a)0(a) . 



(66) 



Expanding eq. (|63j) and comparing with eq. (|H7j) enables an expression to be 
found for of 



for 0(a) improved Wilson fermions and where 

g~ s = go (i + ^iso 2 + • 



(67) 



(68) 



at the scale \i = 1/a. t\ is known and is given by t\ = 0.468201, [10] ■ Hence 
may be computed. Values are given in Table El 

For two loops a simple exact analytic expression is possible for AZq{M) of 



AZ, 



o 



2b g 2 n 



2b 



1 , bl 2 

bo 



b d O*i - b l d O;0 



(69) 



The expression in eq. ()69() is the result of renormalisation-group-improved 
boosted perturbation theory. We can finally tadpole improve it to obtain to this 
order 



AZ^(a) 



Ur 



1—nr 



2b g 



2b 



1 



6i 



gl 



b d O:l -"l d O;0 pi b 
2b b 1 ' 4 6l 



(l-n D ) 



(70) 



12 Appropriate numerical values of Uq = (jrTrt/n) are given in .21 . For other numbers 
a simple interpolation between these values can be performed, or alternatively a Pade fit, 
including the known first three coefficients of the plaquette expansion, can be made. 
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where no is the number of derivatives in the operator. Note that for one derivative 
operators, TI has no effect, ie eq. is the same as eq. (|70p. Again we have the 
factor, eq. (fTHj) . 

Zg*{a) = AZ^(a). (71) 

Thus AZ0 4T (a) is the function that takes you directly from the lattice result to 
the RGI result. 

Finally if we additionally wish to TI the improvement coefficients, [HI], then 
we replace g$ by in eq. (J3*7)l . Numerical results from this procedure for Co 
in v 2 with C2 = are given in Table 0] The associated k c , necessary for the 
computation of am q , eq. (|3l)J) are given in [TFj . 



5.2.3 Non-perturbative Z determinations 

We now look at the RI'— MOM non-perturbative determination of renormalisation 
constants, using the method proposed by Martinelli et al [IT] which mimics (up 
to a point) the approach of the perturbative lattice procedure, by defining 



y RI'-MOM /_\ 

y Rl'-MOMt \ _ q;LAT \P I 

^o-.lat Krpi a ) — ~ 



-Tr 
12 ±L 



Ao(p)^o]born(p) 



(72) 



where the wave function renormalisation constant Z MOM (p) can be fixed from 
the conserved vector current or from the (Fourier transformed) quark propagator 
3 q (p) by 



Tr(-i Ea 7a sm(ap x )aS (] 1 (p)) 
12 Ea sin 2 (ap A ) 



Z q .,Zr- MOM (p) = V 2^Z, q ■ (73) 



(There are still various possibilities for Z g;Mr , see eg jlHj for different definitions. 
Again Ao is the amputated three-point Green's function for the operator O.) For 
our implementation using a 'momentum source' see [49j . For the higher deriva- 
tive operators considered here, this is a non-covariant renormalisation condition, 
depending on the momentum direction, EO] (numerically this is a small effect 
for the momenta considered here). 

The non-perturbative results for Z .^ M0M should now be brought to an RGI 
form, which can only be done perturbatively. In order to avoid problems caused 
by the non-covariance of the renormalisation condition eq. (172)) we first transform 
Zo-lat M ° M (perturbatively) to a covariant scheme S like MS or MOM employing 
a conversion factor of the form 

l + cf(^) 2 + cf (^) 4 + ... • (74) 



For S = MS the general expression to one-loop order can be found in jlU], while 
an explicit formula for V2b is given in jSU].The one- loop expressions for S = MOM 
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are also known. In the case of v 2a and V2b, three-loop expressions can be derived 
from [5T] and will be used in the following. 

In a second step, multiplying the resulting numbers by AZ@, we obtain Zq 31 
(see eq. (|I5|)). Thus AZ@ has to be found, which as in section 15.2.21 is again the 
computation of the perturbative coefficients of the anomalous dimension 7 5 . For 
S = MS the anomalous dimension is known to three loops for v 2 , ^3 and V4. For 
S = MOM only the first two loops (or coefficients) are available in the case of 
t> 3 and t> 4 , while in the case of v 2 we can make use of the three- loop calculation 
in [5T] . For example, expanding in g MS , 7^ = do ]0 (g MS ) 2 + df ) . l (g MS ) 4: + . . ., we 
have to two loops, similarly to eq. (fTUj). 

b d o- 1 - b l d O;0 
26^1 

(75) 

Depending on the choice of S and the coupling in which 7 5 is expanded, there are 
several possibilities (of course equivalent if one knew the whole power series) . We 
shall briefly describe two methods here (I and II), whose difference we shall use to 
estimate the potential error due to (unknown) higher terms in the perturbative 
expansion. In method I we choose S = MS and expand in g MS . In method II 
we work in the MOM scheme (see [M]), and therefore it may be more natural 
to expand in other coupling constants defined using momentum renormalisation 
conditions. In [52] several possibilities are given. We shall use here the coupling 
defined by using the three gluon vertex, the MOMgg scheme, in the notation of 

If we plot Z^ 1 against (r /i p ) 2 and the perturbative expressions are suffi- 
cently well known, we would expect to see a plateau where Z^ 1 is independent 
of (rofj,p) 2 . This region occurs when fi p is not too small ~ A otherwise non- 
perturbative effects play a role, nor when it is too large as lattice artifacts then 
occur. Unfortunately these O(r /i p ) 2 effects may become large (but depend very 
much on the operator considered). In perturbation theory in the chiral limit 0(a) 
terms of the Green's function have opposite chirality to the leading term, they 
disappear when the trace in eq. (|72p is taken. (For explicit lowest order results 
see (Ml-) Condensates may spoil this at low /i 2 , but here we are looking at higher 
momentum scales. Thus we shall take the plateaus in the chiral limit as our 
renormalisation constants. 

We have made determinations of Z .^ M0M using eq. (f72|) for (3 = 6.0, 6.2 
and 6.4 on 24 3 x 48, 24 3 x 48, 32 3 x 40 lattices respectively. For each f3 three or 
more quark masses were used and a linear extrapolation in am q was performed 
to the chiral limit. More details will be given in |53j . 



2b (g MS f 



^1 I _MS\2 



1 + 

bo 
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u.u 


fi 9 


fi 4 




F 1 


1 A 1 P. 

1. 41:10 


1 47^ 

1.4: 1 


1 ^99 


7 RGI 
/j V 2a 


TDD TD'-p 

1 sxD-r 1 


1.536 


1.571 


1.604 




ATT3 T 

INF — 1 


1.46 


1.51 


1.56 




ATT3 TT 

INF — 11 


1 45 


1 5D 


1 55 




11 


1 939 


1 918 


1 90^ 




F 1 


1 4fl7 


1 .4:00 


1.010 


*7 RGI 

Z v 2b 


'I'll 1 1 I v' r ' 

1 KB-F 1 


1.519 


1.555 


1.589 


NP - 1 


1 4<S 


1 .OU 


1 . 00 




INF - 11 


1 4^ 

1 .4:0 


1 4Q 


1 ^4 
1 . 04 


r TI 


TT 

J. J. 


t OAK 


1 99Q 
L.ZZif 


1 91 




PT 


1.928 


2.038 


2.129 


17 RGI 


TRB-PT 


2.268 


2.330 


2.389 


NP - I 


2.2 


2.3 


2.4 




NP - II 


2.1 


2.2 


2.3 




PT 


2.367 


2.548 


2.700 


7 RGI 


TRB-PT 


3.156 


3.242 


3.325 


NP - I 


3.1 


3.3 


3.5 




NP - II 


2.9 


3.1 


3.2 



Table 4: Z mi results (and some cq results) at (3 = 6.0, 6.2 and 6.4, for v n , n = 2, 
3 and 4. PT denotes the perturbative results from section 15.2.11 and TRB-PT from 
section 15.2.21 Note that to obtain Z XI for PT the three-loop results from Table |H1 have 
been used. Results from both NP variations, I and II are shown. 

5.3 Comparison of Z mi results 

We shall now use these results to find Z£* Z™ 1 and Z£» . In Figs. EU [TJJ [TJ 
we plot Z^ , Z^ 1 and Z^f for (3 = 6.2 as computed from these techniques. In 
the case of the non-perturbative Zs a reasonable plateau is seen for large (r /z p ) 2 
enabling a value for the renormalisation constant to be found. Method II seems 
to reach a plateau faster than method I, so we shall use these results (taken 
around fi p = 5 GeV), with the difference to method I giving a rough estimate of 
the error. Also shown in the plots are the TRB-PT results. The NP and TRB-PT 
results lie close to each other, with a maximum discrepancy of about 8%. In the 
NP determination of the Zs, the plateaus become better as j3 increases. This is 
shown in Fig. 

In Table 0] we give the results from PT, TRB-PT and NP method (both 
variants) for various renormalisation constants 13 at (3 = 6.0, 6.2 and 6.4. 

13 For an independent nonperturbative calculation of Z^ 1 see [51) . 
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Figure 10: at = 6.2 versus (r /x p ) 2 . NP method I is denoted by circles, while 
method II is given by filled squares. The dashed line is the TRB-PT result, while the 
full line is the NP estimate from method II. 
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Figure 11: Z^ 1 at f3 = 6.2 versus (ro/j, p ) 2 . The same notation as for Fig. 1101 
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Figure 12: Z^ 1 at = 6.2 versus (ro/^ p ) 2 . The same notation as for Fig. 1101 
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Figure 13: versus (ro/i p ) 2 for (3 = 6.0 (circles), 6.2 (filled squares), 6.4 (diamonds) 
using method II. The corresponding dashed lines are the TRB-PT results, while the 
full lines are the NP estimates. 
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Figure 14: v^ 1 versus (rom ps ) 2 for j3 = 6.0 (circles), 6.2 (squares) and 6.4 (diamonds). 
The filled symbols are obtained when using the TI value for cq while the empty symbols 
have Co = 0. Also shown are linear extrapolations to the chiral limit (dashed lines). 
Other notation as in Fig. 



6 Results: chiral and continuum extrapolations 

6.1 The phenomenological approach 

From the bare results given in Appendix [01 and the Zs in Table 0] we can now 
construct our renormalised matrix elements. As well as the continuum extrapola- 
tion a — > 0, as the quark masses presently used are rather heavy we also need to 
extrapolate these renormalised results to the chiral limit. In this section we shall 
consider both these limits. In Fig. El we show the results for V2b for (3 = 6.0, 

6.2 and 6.4 with both Co = Cq 1 , as given in Table El and Co = 0. (We would 
expect that for Co set to the TI value the results are practically 0(a) improved 
for non-zero quark mass.) 

Also shown in the figure is a linear extrapolation v n = F^(rom ps ) where 

F^\x) =a n x 2 + b n , (76) 
(with n = 2). This fit describes the data well and (not surprisingly) using either 
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Figure 15: v^ 1 versus (rom ps ) 2 for /3 = 6.0 (circles), 6.2 (squares) and 6.4 (diamonds). 
Also shown are linear extrapolations to the chiral limit (dashed lines). Other notation 
as in Fig. H3 



Co value gives the same result in the chiral limit. Indeed the v^b 0(a) improved 
results seem to be independent of the quark mass. 

A similar situation holds for v^b, V2 a , t>3 and V4 all evaluated with non-zero 
momentum but with larger (and increasing) error bars. In Fig. we show the 
results for u 2a for (3 = 6.0, 6.2 and 6.4 together with a linear chiral extrapolation 
(using the TI value for cq). Immediately noticeable when comparing with Fig. fT4l 
is the large increase in the error bars and less consistent ordering of gradients 
with increasing j3. 

In Fig. El we show the results for i> 3 and in Fig. El the analogous results for 
t>4. As expected the results become noisier for increasing (3. Perhaps surprisingly 
the results for V4 seem to be as consistent over our /5-range as those of V3 (we 
have no real explanation for this). 

The last limit to be taken is the continuum limit, a —* 0. As discussed in 
section 13] we believe that for v-i the improvement terms are numerically small and 
so can be neglected. Thus we can make an extrapolation in a 2 (rather than a). 
While we cannot be so confident in this for the higher moments, based on the 
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Figure 16: v^ 31 versus (rom ps ) 2 for (5 = 6.0 (circles), 6.2 (squares) and 6.4 (diamonds). 
Also shown are linear extrapolations to the chiral limit (dashed lines). Other notation 
as in Fig. H3 



experience with the lowest moment, we shall also assume this for these higher 
moments. In Table El we give first the RGI values at f3 — 6.0, 6.2, 6.4, for v^ 1 , 
v^' both for non-zero and zero momentum, v^ 31 and vf 31 . (These have all been 
obtained using the NP method II results, as given in Table EJ) 

We now use these results to perform a continuum extrapolation. In Fig.llSlwe 
plot the continuum extrapolations for the various v-f 3 . A very consistent picture 
is obtained firstly between the different representations ('a' and '6') and secondly 
between the different momenta used in the '6' representation. As expected though 
using a non-zero momentum gives a much noisier signal: in the extreme case 
between the off-diagonal and diagonal representations the error is about ~ 0(2) 
larger. We shall present our final result using v^ 1 for p = only In Fig. 
we show the results for v^ 31 and vf 31 . Using the modified operators of eq. (|3ip 
enables a relatively smooth extrapolation to the continuum limit, giving results 
with about a 20% - 30% error. 

Finally, we convert our results to the MS scheme at a scale of \i = 2 GeV, 
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Figure 17: v^ 31 versus (rom ps ) 2 for (5 = 6.0 (circles), 6.2 (squares) and 6.4 (diamonds). 
Also shown are linear extrapolations to the chiral limit (dashed lines). Other notation 
as in Fig. El 



using the three loop result for [AZjf s (p)] 1 given in Table E3 We find 

wp(2GeV) = 0.245(9), 
vf^{2GeV) = 0.083(17), 

<^(2GeV) = 0.059(18). (77) 

The total error are the combined errors from the three point functions and chi- 
ral/continuum fits together with the error for the renormalisation constant. For 
v 2b this indicates that the dominant error is now possibly coming from the renor- 
malisation constant; the opposite is so for the higher moments. 

How reasonable are the results in comparison with experiment or the MRST 
phenomenological fits? The MRST numbers from section El in Table [0 are also 
plotted in Figs. HH1 and EH We see that for v^ 1 , the discrepancy between the 
experimental result and the lattice result stubbornly remains - and has persisted 
ever since the first pioneering works in the field pQ. It is also notable that v^ 31 
in particular is too large in comparison with the phenomenological result. For 
both vf 31 and v f 31 the chiral and continuum extrapolations are more problematic 
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Figure 18: v^ 1 and v^ 1 (both for p = pi and 0) versus (a/ro) 2 , using the results at 
j3 = 6.0, 6.2, 6.4 (filled circles). A linear continuum extrapolation in a 2 is also given 
(dashed line and empty circle). The star is the MRST value given in Tabled 
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Figure 19: vf 31 and vf 31 versus (a/ro) 2 , using the results at (3 = 6.0, 6.2, 6.4. A 
continuum extrapolation is also given. The same notation as for Fig. 1181 
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p 


P 


6.0 


6.2 


6.4 


oo 


U 2a 


pi 


0.290(37) 


0.379(51) 


0.304(41) 


0.343(56) 


V 2b 


Pi 


0.338(16) 


0.328(16) 


0.340(16) 


0.335(22) 


V 2b 





0.354(8) 


0.344(8) 


0.342(8) 


0.335(11) 


,,H3/ 
y 3 


Pi 


0.141(19) 


0.156(24) 


0.132(20) 


0.137(28) 




Pi 


0.090(14) 


0.099(29) 


0.104(26) 


0.110(33) 



Table 5: v^ 31 results for 0(a) improved fermions at (3 = 6.0, 6.2 and 6.4, for v^ 31 , 
n = 2, 3 and 4. 'oo' denotes the continuum extrapolation (a = 0) limit. 

than for V2b- This can probably only be resolved by increasing the statistic of the 
ensembles and by additional simulations at other (3 values. 

Finally for completeness we give the results for «™ w for q = u, d separately. 
We find 



4" )W (2GeV) 
4 up (2GeV) 

v {u)MS 



(2 GeV) 



0.436(18) 
0.136(25) 
0.096(25) 



(78) 



for the u quark and 



(d)Ms 

(d)MS 



(2 GeV) 
(2 GeV) 



t;f MS (2GeV) 



0.191(7) , 
0.052(11) 
0.027(15) 



(79) 



for the d quark. As discussed previously in section l4~3"l the quark line disconnected 
(qldis) contribution to the matrix element is not computed. Thus on the RHS 
of eqs. (J78|) . (fT9*|) there should be an extra term v^ IS \ q idi S , which is the same for u 
and d quarks, and of course, cancels for the NS results of eq. (J77J). 



6.2 Chiral perturbation theory 

Although linear extrapolations in rrii s seem to describe the results presented 
earlier quite well, it is not clear that for the quark masses used here (and in 
other simulations) higher order effects and/or chiral logarithms can be neglected. 
There has recently been a flurry of interest in this direction. In jH] , based on chiral 
perturbation theory proposed in 0] , a fit function model is used which tries to take 
into account the 'pion cloud' around the nucleon, giving with v^ 31 = F^ n )(r m ps ), 

FW(x) = a n x* + b n (l-cxHn — * V (80) 
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where A x is a parameter, the chiral scale, usually taken to be of 0(47rf n ) ~ 
0(1 GeV) (where f n ~ 93MeV). For A x = or large pseudoscalar mass the 
equation reduces to the linear case, eq. (ff6|) . (Unfortunately as most of the 
masses used in the present simulation are larger than the strange quark mass, 
eq. (fKU|) may need higher order terms in chiral perturbation theory to be valid 
at these large masses.) These first results have been confirmed by further chiral 
perturbation computations, [HJ [7J |H1 E] 14 ■ In particular in quenched QCD was 
considered, with the result that at least for the nucleon there are no additional 
quenched chiral logarithms present, so-called 'hairpin diagrams', so the structure 
of the result in eq. ()8())1 remains unchanged. Furthermore in (unquenched) QCD 
we have 

c =7T^~°- 66 > (81) 

while for quenched QCD, assuming that ~ cr 1 ), (3^ ~ (3^ in [§], then c ~ 
0.28. Of course, in principle these formulae, eq. (J80|) . are valid in the continuum 
so we should first take the continuum limit and then apply chiral perturbation 
theory. Thus we should interpolate the values for v n in Figs. I14 | ITK j fTBl and IT71 to 
a set of constant (rom ps ) 2 for each value of (3. For each of these values of (rom ps ) 2 
a continuum extrapolation should be performed. A chiral extrapolation of the 
data can then be attempted. Unfortunately our 'grid' of data points is not fine 
enough and also for the higher moments is too noisy for this precedure. Thus we 
shall try a 'half way house' approach and attempt a simultaneous continuum and 
chiral extrapolation of the data, 

vf 1 = {r m ps ) + c n 2 + d n ar m 2 ps , (82) 

where the first term represents chiral physics, given by eq. (|7U|) or eq. (jHUjl . the 
second term discretisation effects and the last term accounts for residual 0(a) 
effects oc am q , see eq. ()35|) . With this type of fit the number of free parameters is 
slightly reduced in comparison with the previous fit procedure given in section loTI 
and so tends to produce smaller error bars. We shall restrict our results here to 
our best data set - v^ 1 and first check that using eq. (fTUj) for reproduces 
our previous results. In Fig. EU]we first fit Vzb to eq. (J%2~j) with given by the 
linear function of eq. (f7B|) . F^> = a 2 (r m ps ) 2 + b 2 , and then plot v^ 1 — c 2 (a/r ) 2 — 
d2arorrip S for our three (3 values. (There is no perceptible difference between using 
Co = or Co = 0, for definiteness we show the result for cq = 0.) The points lie 
reasonably well on a straight line, with extrapolated result 0.331(10) consistent 
with our previously obtained value in Table El The alternative possibility of 

14 These works find the leading chiral logarithm behaviour, 

FW(x) ~ & n (l-cx 2 m(z 2 /(r A x ) 2 )), 
which is built into the model of eq. I|80|l . 
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Figure 20: v^ 1 — C2(a/ro) 2 — d 2 arorni s = = a2(rom ps ) 2 + 6 2 (ie using the chiral 
function of eq. (|76jl) versus (rom ps ) 2 , dashed line. Filled circles, squares and diamonds 
represent j3 = 6.0, 6.2 and 6.4 respectively. The empty square represents the chirally 
extrapolated value. Again the MRST phenomenological value of v^ 31 is represented by 
el stcir. 



vf 31 — (F^ (r m ps ) — rf 2 ar mp S ) versus (a/r ) 2 would display the 0(a 2 ) lattice 
discretisation errors. However from Fig. we see that the 0(a 2 ) effects are 
small and this alternative plot just reproduces them again. 

Bolstered by this result, we now try to use F^ from eq. (|50|) as shown in 
Fig. |^] However it is difficult to detect any non-linearities in the data and a 
6-parameter fit (a 2 , b 2 , c, A x , c 2 , d 2 ) fails. We were forced to see if such a 6- 
parameter fit could be plausible, by fixing the chiral limit, b 2 to be the MRST 
phenomenological value given in Tabled (Note that there is no reason that the 
quenched QCD value must be the same as the QCD value, however for many 
hadronic quantities there appears to be little difference between the quenched 
and unquenched QCD values.) As expected, while the fit (dashed) line and the 
numerical results are in reasonable agreement, all of the curvature of the fit takes 
place at small quark mass values where there is no data. Also the fitted parameter 
result for the chiral scale, A x ~ 505(48) MeV, is small in comparison with the 
expected value discussed earlier. For c we found 0.39(6) which lies between the 
unquenched and quenched values. 

So it would seem that any possible non-linearities can only show up at rather 
small quark mass outside the present range of data. Teraflop simulations will be 
necessary to reach more physical pion masses. At present we shall stick to the 
simplest extrapolation possible. 
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Figure 21: — C2(a/ro) 2 — d2arom 2 ps = F% (using the chiral function of eq. (|8l)|)) 
versus (rom ps ) 2 , dashed line. The chiral limit was fixed to the MRST phenomenological 
value. Other notation is the same as in Fig. I2UI 



7 Conclusions 



In this article we have computed v n , n = 2,3, 4. There is still a difference between 
the lattice results and the phenomenological results, particularly apparent for v-i 
where there is a ~ 40% discrepancy. 

We have tried here to narrow down the range of possibilities for the disagree- 
ment, on the experimental side by comparing the global MRST/CTEQ fits with 
the experimental data: there is good agreement. 

On the lattice side, we have 0(a) improved the lowest moment, investigated 
possible mixing operators for the higher moments and discussed renormalisation, 
both perturbatively and non-perturbatively. At least for i> 2 there do not seem 
to be large 0(a 2 ) corrections and so a continuum extrapolation can be reliably 
performed. For the higher moments, we have introduced modified operators for 
t>3, t>4, which reduces the numerical noise and improves the signal. It is difficult to 
see if there are any 0(a) corrections. However at present the simplest assumption 
that these corrections are small also fits with the data. 

Although this is only a partial study here of mixing operators for the higher 
operators (due in part to the present difficulty of even defining renormalisation 
constants for several of these operators) presently we find little sign of problems. 
Indeed even for v™ 3 the mixing with a lower dimensional operator appears harm- 
less. In other situations this is not the case and there may be significant changes, 
see eg |5T] . 

We have discussed and compared various renormalisation procedures rang- 
ing from simple perturbation theory to TRB perturbation theory to a non- 
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perturbative method. Using NP results eliminates uncertainty concerning the 
renormalisation constant. It is seen though that between TRB perturbation the- 
ory and non-perturbative results there can be up to an 8% difference - far less 
than the difference between the lattice and phenomenological result for v^ 31 . 

Finally although at present we see little numerical evidence of chiral loga- 
rithms, this is perhaps telling us that we must go to significantly smaller quark 
masses before the chiral extrapolation 'bends' over. However preliminary re- 
sults at lighter quark mass for unimproved Wilson fermions also show a linear 
behaviour and so do not seem to improve the situation, jHEl EE] • 

Whether quenching effects are significant remains unclear, but recent un- 
quenched results, [HB1EH!, do not seem to reveal any significant differences be- 
tween quenching and unquenching, at least in the quark mass range considered. 
Finally there are hints of a different situation for overlap fermions, [SU], which 
might suggest again that one has to simulate at light quark masses close to the 
chiral limit - a challenge for the lattice. 
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Appendices 

A Some values of [AZJf*(/i)]~ 





one-loop 


two-loop three-loop 




[A^f(zi)]- 1 


2.00 GeV 

2.12 GeV (1/a at /? = 6.0) 
2.90 GeV (1/a at /3 = 6.2) 
3.85 GeV (1/a at (3 = 6.4) 


0.783(10) 
0.776(10) 
0.743(8) 
0.718(7) 


0.721(8) 
0.715(8) 
0.688(7) 
0.668(6) 


0.732(9) 
0.726(8) 
0.696(7) 
0.674(6) 




[A^f(Zi)]- 1 


2.00 GeV 

2.12 GeV (1/a at (3 = 6.0) 
2.90 GeV (1/a at /3 = 6.2) 
3.85 GeV (1/a at (3 = 6.4) 


0.682(14) 
0.673(13) 
0.629(11) 
0.596(9) 


0.596(10) 
0.589(10) 
0.555(8) 
0.529(7) 


0.609(11) 
0.602(10) 
0.564(8) 
0.537(7) 




[a^ s (a*)] _1 


2.00 GeV 

2.12 GeV (1/a at (3 = 6.0) 
2.90 GeV (1/a at (3 = 6.2) 
3.85 GeV (1/a at (3 = 6.4) 


0.619(15) 
0.609(15) 
0.559(12) 
0.522(10) 


0.520(11) 
0.512(10) 
0.476(9) 
0.448(8) 


0.534(13) 
0.526(12) 
0.486(10) 
0.456(8) 






2.00 GeV 

2.12 GeV (1/a at /3 = 6.0) 
2.90 GeV (1/a at (3 = 6.2) 
3.85 GeV (1/a at (3 = 6.4) 


0.268(10) 
0.261(9) 
0.228(7) 
0.205(6) 


0.195(6) 
0.191(5) 
0.170(5) 
0.156(3) 


0.201(6) 
0.196(6) 
0.174(5) 
0.159(4) 



Table 6: Useful values of [AZ£ fl (fi)}- 1 (n = 2, 3, 4) and of s {n) = {g MS \n)) 2 /4tt. 
The errors are a reflection of the error in A MS ro. The lattice simulations performed 
here, as described in section give the above 1/a values (found from using r^/a in 
38 , namely r^/a = 5.37, 7.36, 9.76 at (3 = 6.0, 6.2, 6.4 respectively, together with the 
scale ro = 0.5 fm). 
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B Group properties of mixing operators - an 
example 

When QCD is put on the lattice, the necessary analytic continuation from Mink- 
owski to Euclidean space replaces the Lorentz group by the orthogonal group 
0(4), which by the discretisation of space-time is further reduced to the hyper- 
cubic group H{4) C 0(4). Since H{4) is only a finite group, the restrictions 
imposed by symmetry are less stringent than in the continuum and the possi- 
bilities for mixing increase, sometimes in a way which is not easily anticipated 
by our (continuum) intuition. In this appendix we illustrate the problem by an 
example. For a more complete and systematic treatment we refer to 

As a further symmetry we have charge conjugation. It operates on the fermion 
fields q(x), q(x) and on the lattice gauge field U^{x) according to 

q(x) 2* Cq(x) T , 
q(x) Z -q{x) T C~\ 
Up{x) $ U^x)* (83) 

with the charge conjugation matrix C satisfying 

CllC- 1 = - 7 , • (84) 



So we get, eg, 



(9 7 (— 1 ) n (1' y 



Identifying elements R of H(4) with the corresponding 4x4 matrices in the 
defining representation we find that R acts on qj^q as follows: 

Ti^ Rv^qiuq , (86) 



V 



ie, the four operators q^y^q form a basis for the defining representation of H(4). 
More generally, we get for the action of R: 

Ul,U2,...,U n 

and 

(', 

fllfl2---fJ. 



f0775 _R A p f D D D . . . D (H775 



I/l,i/ 2 ,...,l*n 



Thus the 4 n operators 0^ lft2 _ „ as well as the operators C77L...ii„ f° rm a basis for 
a representation of H(4), which for n > 1 is reducible. It is helpful to consider 
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these operators as forming an orthonormal basis of the representation space. 
(Orthonormal) bases of irreducible subspaces have been given in [35]. 

Operators transforming according to the same irreducible representation of 
H(A) (and having the same C-parity) will in general mix with each other so 
that one has to consider appropriate linear combinations. Writing down these 
linear combinations one has to choose the bases in the different (equivalent) 
representation spaces such that they transform identically under H(A) and not 
just equivalent ly. 

Consider two bases which are known to transform according to the same irre- 
ducible representation. In order to check whether they even transform identically 
it is sufficient to compare their transformation behaviour under a set of genera- 
tors. For H(A) there is a set of three generators {a, (3, 7} given by 

/ 1 ^ 
10 

10 

V 1 J 

( 1 ^ 

1 
1 

V 1 



a 



7 



/ -1 



V 




1 





1 



J 

°\ 



1 / 



According to (J87|) . a interchanges 1 and 2, j3 performs a cyclic permutation of 
the values of the indices, and 7 produces a factor of —1 for every index 1, eg, for 
the operator Oj ul one finds 



CV 

^3341 


a 


CV 

^3342 > 


CV 

^3341 


jr. 


^4412 > 


^3341 


7 


-CV 

^3341 



(90) 

As an example for mixing operators, let us consider VA . It belongs to a 
doublet of operators transforming according to the H(A) representation (in 
the notation of [55] and Table |2J) and has positive C-parity. Indeed, the two 
operators 



y^. (cyy _l (-97 

9 V {1122} + ^{3344} 



CV 

^{1133} 



CV 

^{2244} 



and 



«2 



1 

71 



ocv + ?o y — cv — n 

ZLy {1144} + ZLV {2233} ^{1122} U 



_ CV - CV 

{3344} ^{1133} ^{2244} 



(91) 

V20 V4 
(92) 
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form an orthonormal basis for this representation. It is obvious that both oper- 
ators remain unchanged under the action of 7. The generator a acts on them as 
follows: 



U 2 — > 

while (5 gives 

ui \ui + ^V3u 2 , 
P 1 

U2 -> 3 



-|>/3ui - , (93) 



/3 



^ui-^- (94) 



lUi := 



It is straightforward to check that the operators 

( ~~ ^1122 ~~ ^2112 ~~ ^1221 ~~ ^2211 + 2Ci212 + 202121 

_ r> 7 — n 7 — n 7 — n 7 -i- 9n 7 -i- 9n 7 

^3344 ^4334 ^3443 ^4433 + ZLy 3434 + ZLy 4343 
+ ^1133 + ^3113 + ^1331 + ^3311 — 2O1313 — 2C?3 131 
+ C^2244 + O4224 + ^2442 + ^4422 — 202424 — 2(^4242) , (95) 



and 



w 2 — ( + ^1122 + ^2112 + ^1221 + ^2211 201212 202121 

+ 03344 + 04334 + 03443 + 04433 — 2O3434 — 2(94343 

+ 01133 + 03113 + 01331 + 03311 — 2(9i3i3 — 20 3131 
_l_ CV 4- CV 4- CV 4- CV — 9H 7 — 9H 7 

1 ^2244 ' ^4224 ' ^2442 ' ^4422 z,L/ 2424 z,L/ 4242 

— 2C>7i44 — 202 114 — 2O1441 — 2(92 411 + 4(97414 + 404 141 

~~ 202233 ~~ 2O3223 — 202332 ~~ 2O3322 + 402323 + 40323 2 ) , 

= \0Z\ (96) 
transform identically to Ui, u 2 : 

w 1 A Iiyj - \\fZw 2 , 



U>2 



-|a/3wi- |w 2 , (97) 



u> 2 A iy/3 Wl -±w 2 , (98) 

and they are invariant under the action of 7. Hence any linear combination 
r ■ Ui + s ■ Wi (i = 1,2) has the same transformation properties under H(A) as 
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U{. In particular, V4 may mix with O™ 1 , and the renormalisation of V4 will in 
general involve also O™ 1 . 
Also the operators 



1 

^ ( + ^1324 — ^2314 — ^1423 T" ^2413 T" ^3142 — ^4132 ~~ ^3241 "I" ^4231 

I /nTT8 /o775 //-)775 , /n778 , /nTYB /n77s /nTTs j_ /n775 ^ ('QQ'l 
T L/ 1234 — ^3214 — ^1432 "I" ^3412 ~+~ ^2143 — ^4123 — ^2341 ~+~ ^4321 j i l y J ) 



„,,' / I /n775 /nTTS /oTTs _L /*Yns i /n775 /»1776 /n775 , /n775 

•— . I ~T ^1394 ^9314 ^1493 ~r ^9413 I ^3149 ^4139 ^3941 ~r <S . 



and 



' ^ f i p\nh _ /nrrs /nrys i /otts _l ^otts /httb /otts _l /n776 

_ 4v^3 34 3214 ^1432 i* ^3412 i* ^2143 ^4123 ^2341 i* ^4321 

/*YY75 _L POTT'S _L /*Yns /nTTS _L /»1TY5 _i_ /n775 /nTTB 

^1324 ~r ^9314 ~r ^1493 ^9413 ^3149 ~r ^4139 ~r ^3941 



1324 ~ ^2314 ~ ^1423 ^2413 ^3142 ~ ^4132 ~ ^3241 ^4231 
-|- ^L^i243 ZLy 4213 ZLy 1342 "t" ZLy 4312 i* ZLy 2134 ZLy 3124 

ZLy 2431 + ZLy 3421, 

' 0™\ (100) 



2V3 



transform in exactly the same way under a, (3, and 7, ie under H(4), as follows 
from (|88p. Therefore (9™ 2 may mix with G U4 , too. 

C Two- and three-point correlation functions 
C.l General formulae 

In this appendix we shall give explicit expressions for the correlation functions in 
terms of quark propagators. We start with the two-point correlation function. A 
suitable proton operator is 

B a (t; p)=Y< e- l ^e abc u a a {x, t) [u b {x, t) TD C-f 5 d c {x, £)] , (101) 

X 

(a . . . denote colour indices, a . . . Dirac indices, T D is the transpose in Dirac space 
and 75 = 71727374)- C, the charge conjugation matrix, has the defining property 
given in eq. (|83jl and is antisymmetric. As we take our gamma matrices to be 
hermitian then C may be taken as as unitary. Thus — C = C T ° and = . 
One possible choice (used in our computer programme) is C = 7472 (so that 
— C = C^ 1 ) and the Dirac basis 

(but the results given here should not depend on this particular choice). The 
last two quark fields in eq. (jlOlj) form a di-quark structure, while the first quark 
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carries the spin index. With the help of 



B a {t;p) =Y / e l ^e abc lf a (x,t) [d\x, t)-y 5 Cv°(£, t) 

X 

the two-point correlation function is formed in the usual way 

C r (t;p)=Ti D T(B(t,p)mp)), 



(103) 



(104) 



where we have introduced an arbitrary Dirac matrix, T, which for unpolarised 
nucleons, the case considered here, should be taken as |(1 + 74). This equation 
may be re-written using quark propagators propagating from a (point) source to 
a sink. As we are averaging over the gauge fields we may shift all sources to (0, 0). 
Some algebra then yields 

C T {t;p) = -V s Y,e- i ^e abc e a ' h ' c ' 



Tr 



D 



TG^ aa '(x, t; 0, 0)1 Tr D \G {d)bb ' {x, t; 0, 0)G^ cc '{x, t; 0, 0) 



+ Tr 



D 



TG {u)aa> t . Q ; ty G (d)bb> ( - t . Q; ty G {u)cc> ( - f . ^ q) 

where we have defined a tilde in Dirac space by 

X = {C l5 X l5 C) T ° . 



{U} 



(105) 



(106) 



The problem is thus reduced to finding the propagator or Green's function for 
quark q from a source (0, 0) to (x, t). In general, the quark propagator from y to 
x is defined by 

and can be computed from 



(107) 



J2M^{x;w)G {q) {w;y) = S : 



xy ; 



'108) 



where M®a is the Wilson (clover) matrix, given in Appendix O 
For the three-point correlation function, 

C r (t,r;p; O q ) = Tr D T(B(t,p)O q (r)B(0,p)) , 



(109) 



we shall only consider the quark line connected term (ie the left diagram of Fig. HJ . 
First we re-write the operator insertion generally as 



O q {r) 



y 

Y,t a {v)Otp{v,w ] y,r)q%w) 



(110) 



y,v,w 
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(ie sum over spatial planes, where (y,r) is the 'center of mass of the operator'). 
For operators without derivatives and with exactly one derivative we have q^yq 
and ^q^^D — D)q for O q respectively, (D and D are defined in eq. (jSOJ)) while for 
two and three derivative operators, to minimise the extension on the lattice we 
'integrate by parts' and choose —qj DD q and \q"i(DDD — DDD)q respectively 
This also allows the higher derivative operators to be built from the previously 
constructed lower derivative operators. 

Some algebra yields the result for eq. ()109j) of 

C r (t,r;p;O q ) 
= -v s E ( Tr ^ 

y,v,w 

where s|? ^(0, 0; v; p, t) is given by 

4 9) (0, 0; v; p,t) = J2 4 9) & t; 0, 0; p)G<«> (x, t; v) , (112) 

X 

in terms of 

Sp )afa (x, t; 0, 0;p) = e - i P-* e abc e a ' b ' c ' x (113) 
G Wbb '{x, t; 0, 0)G {u)cc '(x, t; 0, 0)r + Ti D [G {d)bb ' (x , t; 0, 0)G (u)cc '(f , t; 0, 0)]r+ 

TG (u)bb '{x, t- 0, 0)G (d)cc '(f, t; 0, 0) + Ti D [TG^ )bb '(x, t; 0, 0)}G (d)cc ' (x , t; 0, 0)]] , 

when q = u and a slightly simpler expression for namely 
S^ d)a ' a {x, t- 0, 0;p) = e -^ e ^c e a'b'c' x ( U4 ) 
~& u)bb ' (x, t; 0, 0)f G (u)cc ' (x, t; 0, 0) + Tr D [rG (u)f *' (x, 0, 0)G {u)cc ' (x, t; 0, 0)]" . 

Practically we must find Ej? ^ from a second Green's function using these rather 
ugly looking expressions as sources. By considering 

E 4 9) (0, 0; v; p, t)M® (v; v') = (v't; 0, 0; 0)6 <t , (115) 

V 

we see that this is the wrong way around for the inversion in eq. (jlU8j) but taking 
colour/spin components and using M^(x;y) = (7 5 M*(?/; x) ba ^^)p a the equation 
for E can be re-written as 

E (t/; v^f™ (0, 0; v- p, t) = 75 sf (v f , t; 0, 0; p)6 v ,, t , (116) 

V 

where \cd is the Hermitian conjugate in colour and spin space. We see that in 
this form S is rather like a Green's function from the source given on the RHS 
of the equation. 



4 9) (0, 0; v; p, t)0(v, w; y, r)G® (w; 0, 0) 



an 



111) 
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Thus finding the three point correlation functions is a two step process: first 
the usual Green's function from (0, 0) to any point x is found, stage I. Then a 
second inversion (stage II) is made with the source given from eq. (jll3)) if the 
inserted operator consists of u quarks or using eq. (jll4j) for d quarks. 

The advantage of this procedure is that by tying together the two Green's 
functions appropriately any operator can be inserted with no additional com- 
putational cost. (There is no restriction on the derivative structure and Dirac 
matrix 7). The disadvantage is that for each nucleon state (ie V determining 
whether the nucleon is unpolarised or polarised) , momentum p and nucleon sink 
position (t) a separate inversion is required. Thus results for a range of t values 
are expensive and practically we must restrict ourselves to one value. 



C.2 The Non-Relativist ic Projection 

To improve the overlap with the nucleon we have used Jacobi smearing, [31] , and 
non-relativistic, NR, projection of the nucleon wavefunction. For completeness 
we now briefly describe our implementation of this projection. Rather than using 
the nucleon operator of eq. (jl01|) we shall replace it by 



X 



-ip-x ^abc^ ( a 



u a a (x,t) 



{x,t) (C 75 i(l + 74))^^(f,t) 



(117) 



ie we replace the matrix C75 — > C^§\{\ + 74). Both operators (eqs. ()117|) and 
(|101j) ) behave the same way under rotations and reflections in the spatial direc- 
tions, and both have the same quantum numbers (colour neutral, spin |, isospin 
\ and parity +), and both will therefore overlap with the proton. 

As we shall see, not only is B NR computationally cheaper by a factor two, but 
it has a better overlap with the proton. This can be easiest shown if we use the 
Dirac basis eq. ()102|) first giving 



C75 = 747275 



/ 1 ^ 

-10 

1 

V 0-10/ 



Ct4(1 + 7 4 ) 



( 10 0^ 
-10 


V 000/ 



If we now write out eq. (jll7|) for a spin-up (ie a = 1) proton we have 

B m = t abc (^a u b d c _ u a u b d c 



'118) 



(119) 



setting p = and suppressing the co-ordinate index for simplicity. This is not 
quite the final form. When we sum over all fermion line diagrams, only the 
part of the operator which is antisymmetric under the interchange of the two 
u quarks makes any contribution to the measured Green's function. Since the 
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colour wavefunction e abc is completely antisymmetric, this means that the part 
of eq. ()119|) which survives is 

B NR = e abc (^a u b d c _ lyayb^ _ ly^yb^ . ( 12 0) 

Since the mid-sixties it has been known that the lowest-lying octet and decuplet 
baryons are well described as a 56 of SU (6). According to 577(6) the flavour/spin 
wavefunction of the spin up proton is (eg |61| ) 



V] = V I {2ufUfdi - u T ujd T - M|M T rf T ) . (121) 

Comparing this with eq. ()120|) we see that this is exactly the wavefunction we have 
been using. The wavefunction of eq. (|121|) makes some very successful predictions, 
for example that the ratio of proton to neutron magnetic moments should be — |. 
The experimental value is —1.460, in good agreement with the SU(6) prediction. 
This success suggests that eq. (|120|) is close to the true proton wavefunction. We 
would expect this wavefunction to work even better with heavier quarks, so it 
is an appropriate wavefunction to use on the lattice. If we carry out the same 
exercise with eq. (110 lj) we obtain 



+\u a x u\dl + \ulu\d\ - \u\u\dl - \utu\dl) , (122) 



B x = e ab \u\u\d c 2 - \u\u\d\ - \ulu\d\ 

1 ii a ii b rl c -I- 1 ii a ii b rl c 

^13 4 ' ^314 2"'1"'4"'3 2 

as the explicit component expression. This has many terms involving the 3rd 
and 4th Dirac components with amplitudes just as large as the terms with upper 
components only. In any sort of constituent quark model we would expect these 
terms to be small in the ground state. Adding them into the wavefunction not 
only increases the cost of the computation, it also degrades the signal by adding 
terms which are likely to have more overlap with excited baryon states. 
The correlation functions are constructed also using 

X R (t;p) = J2^ ahc <&t) [3 6 (x,t)i(l + 74)75^^,^] , (123) 

X 

ie replacing 75C — > |(1 + 74)75(7. Thus the NR projection can be obtained by 
projecting out the positive eigenvalues of 74, ie by replacing each quark field, q, 
by 

g^i(i + 74)g, g^gi(i + 74 ) ? (124) 
everywhere and considering polarisation matrices T which satisfy 

i(i + 74)r = r = ri(i + 7 4). (125) 

This gives the NR nucleon two-point function, eq. (|104j) . (So for the Dirac basis, 
eq. ()102|) . only the components a = 1, 2 are non-zero, as discussed above.) In 
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-io 

5 10 15 20 25 30 35 40 45 

t/a 

Figure 22: Relativistic (open circles) and non-relativistic (filled circles) two-point 
operators, eq. (|104j) with p = 0, for (3 = 6.2 at re = 0.1344. 

Fig. 1221 we show a comparison for the nucleon two-point correlation function 
using the relativistic and non-relativistic operators when applied to eq. (jl()4j) . 
The gradient of the left branch measures the nucleon mass. It can be seen that 
this branch has been extended by about 5 units of t/a when using the NR operator 
as opposed to the relativistic operator. 

For the three-point functions the new tilde, replacing eq. ([106)1 . is given by 
X = (C 75 i(l + 74)^(1 + 74)7 5 C) Td , (126) 

and obeys 

i(l+7 4 )X = X = Xi(l + 7 4 ). (127) 
Considering eqs. ffTT3|). (fTT4"j) then eqs. (IT25I) . (IT27j) imply that 

|(1 + 74 )S (,?) = = S (9) ±(l + 74) • (128) 

This identity reduces the number of independent components from 4 to 2 for the 
source for the stage II inversion, eq. ()116j) . Again, for the Dirac representation 
eq. (jl()2|) only the first two components are needed. 
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D Tables 



The action used here (in the quenched limit) is 

S=i/3^Tr c Re[l-f/°]+a 4 £ q(x)M^(x; y)q(y) , (129) 

n xy;q=u,d 

where U° is the product of links around an elementary plaquette in the fj, — v 
plane and the Wilson (clover) fermion matrix is given by 

Y / q(x)M^(x;y)q(y) = 

xy 

r 1 K 

J2{-q( x )(i( x ) — E v( x ) u l( x - a P) I 1 + 7aJ q( x - a P) 

x a a ^ 

J2q( x ) U h( x ) I 1 - 7m] q( x + ap) 

a fi 

-2nac sw g E \q{x) ^F^ over (x)g(x)} , (130) 

where the hopping parameter, k, is related to the (bare) quark mass via eq. (j3Hj) . 
and we are taking mass degenerate u and d quarks. In eq. ()130|) the quark fields 
are normalised according to the lattice conventions ie they correspond to the con- 
tinuum fields by rescaling q — > \j\plkq. (This introduces a further factor 1/(2k) 
on the RHS of eq. (J3TJJ) when using the raw output for the two- and three- point 
functions.) The last term in eq. (jl3())l . sufficient for on-shell 0(a) improvement 
(with a to be determined function c sw (gQ)) has a 'clover' field strength tensor 
given by 

Ftr\ x ) = E [W - > ( 131 ) 

8ig a 2 ±^ v v * J 

where we have extended the definition of the plaquette, so that the /z, v directions 
can be negative. 

Thus we can re-write eq. ()13()j) as 

J2q( x )M^(x;y)q(y)^ 

xy 

E{^»($ +m q )q(x) + m 0c q(x)q(x) - \a^l{ x ) A ^ A t l( x ) 

X fl 

c sw g E \t x ) cr^r r ( x M x )} , (132) 



4 



where 



A+g(x) = -[U li (x)q{x + ajl)-q(x)\ 



» - a 



A q(x) = - \q(x) - UUx - ajl)q(x - ap,)] , (133) 
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K 


Volume 


# configs. 


am ps 




6.0 


1.769 


0.1320 


n 

16 3 x 32 


0(445) 


0.5412(9) 


0.9735(40) 


6.0 


1.769 


0.1324 


16 3 x 32 


0(560) 


0.5042(7) 


0.9353(25) 


6.0 


1.769 


0.1333 


16 3 x 32 


0(560) 


0.4122(9) 


0.8241(34) 


6.0 


1.769 


0.1338 


16 3 x 32 


0(520) 


0.3549(12) 


0.7400(85) 


6.0 


1.769 


0.1342 


16 3 x 32 


0(735) 


0.3012(10) 


r\ r-7 r\r\ / ,4 o \ 

0.7096(48) 


6.2 


1.614 


0.1333 


24 3 x 48 


0(300) 


0.4136(6) 


0.7374(21) 


6.2 


1 CIA 

1.614 


0.1339 


O A 3 * , A O 

24 X 48 


C(300) 


0.3565(7) 


0.6655(28) 


6.2 


1.614 


0.1344 


24 3 x 48 


0(300) 


0.3034(6) 


0.5963(29) 


6.2 


1.614 


0.1349 


24 3 x 48 


0(470) 


0.2431(6) 


0.5241(39) 


6.4 


1.526 


0.1338 


32 3 x 48 


0(220) 


0.3213(8) 


0.5718(28) 


6.4 


1.526 


0.1342 


32 3 x 48 


0(120) 


0.2836(9) 


0.5266(31) 


6.4 


1.526 


0.1346 


32 3 x 48 


0(220) 


0.2402(8) 


0.4680(37) 


6.4 


1.526 


0.1350 


32 3 x 48 


0(320) 


0.1923(9) 


0.4156(34) 


6.4 


1.526 


0.1353 


32 3 x 64 


0(260) 


0.1507(8) 


0.3580(47) 



Table 7: Parameter values used in the simulations, together with the measured pseu- 
doscalar and nucleon masses. Note that the statistics refers to the number of (indepen- 
dent) configurations used for the 3-point functions; the masses have sometimes been 
determined with a larger statistic. 

so that (see eq. (jHDJl) 

D,= I (A+ + A;) . (134) 

am q is defined in eq. (JBTj)) and 

am 0c = 1 (1 - 8) . (135) 

In this latter form eq. (|132j) shows the additional 0(a) operators most clearly. 

In Table 13 we give our parameter values used in the (quenched) fermion sim- 
ulations together with the pseudoscalar, am ps , and nucleon mass, amjv. 

We now give a series of tables tabulating the bare matrix elements for 
q — it, d, the improvement operator matrix elements av^'^ for i — 1, 2 and the 
mixing operators v^' for i — 1, 2 and t^ 9 '™ 1 '' for i — 1, 2 and v^'" 13 ^ /a. 
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K 


0.1320 


0.1324 


0.1333 


0.1338 


0.1342 


P = Pi 


°2a 
V2a 


0.359(19) 
0.1696(89) 
0.189(11) 


0.348(16) 
0.1557(82) 
0.1925(95) 


0.335(21) 
0.144(11) 
0.192(13) 


0.382(45) 
0.163(19) 
0.219(31) 


0.312(35) 
0.123(17) 
0.188(25) 


(u;i) 
av 2a 

(d;l) 
av 2a 

(1) 
av 2a 


-0.0254(17) 
-0.01321(93) 
-0.01211(98) 


-0.0261(15) 
-0.01372(85) 
-0.01238(91) 


-0.0288(40) 
-0.0136(12) 
-0.0117(14) 


-0.0288(40) 
-0.0115(21) 
-0.0170(29) 


-0.0226(36) 
-0.0130(25) 
-0.00974(292) 


(u;2) 
av 2a 

(d;2) 
av 2a 

(2) 
av 2a 


0.0190(14) 
0.01025(82) 
0.00872(79) 


0.0162(13) 
0.00933(75) 
0.00693(82) 


0.00676(165) 
0.00592(104) 
0.00111(117) 


0.00188(206) 
-0.00017(142) 
0.00204(200) 


-0.00361(302) 
0.00136(227) 
-0.00448(260) 


p = 


(u) 

V 2b 

(d) 

V 2b 
V2b 


0.4066(32) 

0.1894(17) 
0.2171(18) 


0.4108(27) 
0.1886(16) 
0.2221(16) 


0.4130(48) 
0.1851(26) 
0.2278(29) 


0.4196(78) 
0.1845(41) 
0.2350(52) 


0.414(11) 

0.1798(55) 
0.2336(69) 


(d;l) 
av 2b 

(1) 
aV 2b 


-0.02844(39) 
-0.01490(26) 
-0.01354(21) 


-0.02895(40) 
-0.01503(27) 
-0.01390(25) 


-0.03002(81) 
-0.01585(51) 
-0.01406(54) 


-0.0328(13) 
-0.0185(11) 
-0.0141(11) 


-0.0311(25) 
-0.0189(16) 
-0.0117(16) 


(u;2) 
aV 2b 

(d;2) 
aV 2b 

(2) 
av 2b 


0.03954(52) 
0.02015(31) 
0.01936(29) 


0.03552(52) 
0.01822(35) 
0.01722(32) 


0.0266(11) 
0.01494(72) 
0.01135(69) 


0.0235(19) 
0.0173(43) 
0.00741(162) 


0.0173(43) 
0.0152(30) 
0.00108(219) 


P = Pi 


(u) 
V 2b 

(d) 
V 2b 

V2b 


0.3918(54) 

0.1842(30) 
0.2076(31) 


0.3959(49) 

0.1807(28) 
0.2152(32) 


0.3908(88) 
0.1727(49) 
0.2180(59) 


0.417(24) 

0.186(13) 
0.231(13) 


0.370(19) 

0.155(11) 
0.215(14) 


(u;l) 
av 2b 

(d;l) 
aV 2b 

(1) 
aV 2b 


-0.02642(60) 
-0.01420(37) 
-0.01222(42) 


-0.01472(41) 
-0.01488(42) 
-0.01320(39) 


-0.0288(12) 
-0.01563(85) 
-0.01297(90) 


-0.0266(32) 
-0.0171(20) 
-0.0102(24) 


-0.0268(41) 
-0.0185(47) 
-0.00898(402) 


(u;2) 
av 2b 

(d;2) 
av 2b 

(2) 
aV 2b 


0.03714(80) 
0.01944(48) 
0.01777(52) 


0.03357(78) 
0.01944(48) 
0.01604(50) 


0.0244(17) 
0.01433(11) 
0.00988(109) 


0.0168(45) 
0.0149(26) 
0.00353(335) 


0.0114(63) 
0.0168(58) 
-0.00315(519) 



Table 8: The bare results for V2 from eq. (|3U)) for (3 = 6.0, c sw = 1.769. 
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K 


0.1333 


0.1339 


0.1344 


0.1349 


P = Pi 


(u) 

v 2a 
V2a 


0.403(24) 
0.187(11) 
0.217(14) 


0.393(23) 
0.179(11) 
0.214(13) 


0.413(38) 
0.182(17) 
0.232(23) 


0.422(42) 

0.175(18) 
0.246(28) 


(u;l) 
av 2a 

(d;l) 
av 2a 

(1) 
av 2a 


-0.0252(18) 
-0.01294(96) 
-0.0122(10) 


-0.0247(18) 
-0.01224(97) 
-0.0125(11) 


-0.0251(28) 
-0.0130(16) 
-0.0122(19) 


-0.0269(32) 
-0.0123(18) 
-0.0144(22) 


(u;2) 
av 2a 

(d;2) 
av 2a 

(2) 
aV 2a 


0.0117(12) 
0.00622(71) 
0.00550(71) 


0.00453(116) 
0.00308(68) 
0.00143(81) 


-0.00110(179) 
-0.0067(110) 
-0.00140(146) 


-0.00493(254) 
-0.00224(144) 
-0.00273(209) 


p = 


(u) 
V 2b 
(d) 

V 2b 
V2b 


0.4100(28) 
0.1920(14) 
0.2179(17) 


0.4020(39) 
0.1839(21) 
0.2181(22) 


0.4071(52) 

0.1834(28) 
0.2237(34) 


0.4052(58) 
0.1775(35) 
0.2278(39) 


(«;1) 
av 2b 
(d;l) 

(1) 
« U 2fe 


-0.02554(35) 
-0.01316(21) 
-0.01239(20) 


-0.02494(48) 
-0.01302(31) 
-0.01191(28) 


-0.02775(65) 
-0.01444(42) 
-0.01323(42) 


-0.02883(96) 
-0.01523(68) 
-0.01357(68) 


(m;2) 
« U 2fe 

(d;2) 
« W 2fe 

(2) 
aU 2fe 


0.02801(41) 
0.01443(26) 
0.01360(24) 


0.02091(58) 
0.01141(38) 
0.00949(35) 


0.01874(79) 
0.01071(56) 
0.00777(58) 


0.0152(12) 
0.00959(87) 
0.00548(97) 


P = Pi 


(u) 
V 2b 
(d) 

V 2b 
V2b 


0.4153(56) 
0.1966(29) 
0.2186(33) 


0.3987(58) 
0.1840(32) 
0.2146(40) 


0.419(14) 
0.1938(77) 
0.2255(81) 


0.406(19) 
0.186(11) 
0.220(10) 


(u;l) 
av 2b 

(d-1) 
aV 2b 

(1) 
aV 2b 


-0.02536(49) 
-0.01316(30) 
-0.01219(30) 


-0.02377(76) 
-0.01251(59) 
-0.01144(42) 


-0.0269(11) 
-0.01438(63) 
-0.01233(82) 


-0.0288(19) 
-0.0148(14) 
-0.0136(15) 


(m;2) 
aV 2b 

(d;2) 
« W 2fe 

(2) 
OW 2fe 


0.02761(67) 
0.01437(42) 
0.01326(44) 


0.01973(92) 
0.01088(74) 
0.00901(56) 


0.0174(15) 
0.0105(10) 
0.00664(125) 


0.0151(25) 
0.00970(179) 
0.00454(215) 



Table 9: The bare results for V2 from eq. (|3U)) for f3 = 6.2, c sw = 1.614. 
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K 


0.1338 


0.1342 


0.1346 


0.1350 


0.1353 


P = Pi 


°2a 
V2a 


0.368(21) 
0.173(11) 
0.196(12) 


0.348(30) 
0.161(15) 
0.185(17) 


0.340(28) 
0.150(12) 
0.191(18) 


0.339(34) 
0.140(16) 
0.199(23) 


0.329(63) 
0.122(24) 
0.205(44) 


(u;i) 
av 2a 

(d;l) 
au 2a 

(1) 
au 2a 


-0.0210(13) 
-0.01092(75) 
-0.01010(73) 


-0.0181(17) 
-0.0092(10) 
-0.00853(99) 


-0.0208(18) 
-0.0104(11) 
-0.0104(11) 


-0.0214(24) 
-0.0103(14) 
-0.0107(16) 


-0.0242(55) 
-0.00981(291) 
-0.0141(39) 


(«;2) 
av 2a 
(d;2) 

(2) 
au 2a 


0.00470(68) 
0.00323(45) 
0.00147(44) 


-0.00203(101) 
0.000043(690) 
-0.00206(64) 


-0.00181(99) 
0.000384(646) 
-0.00229(75) 


-0.00525(152) 
-0.000831(942) 
-0.00472(121) 


-0.00456(307) 
-0.00049(205) 
-0.00352(269) 


p = 


(d) 
W 2b 

V2b 


0.4088(27) 
0.1914(14) 
0.2174(17) 


0.4140(47) 

0.1881(23) 
0.2259(32) 


0.4017(47) 
0.1842(24) 
0.2174(31) 


0.3951(70) 
0.1781(39) 
0.2169(47) 


0.397(13) 
0.1674(71) 
0.2294(90) 


OV 26 

(d;l) 
au 26 

(1) 
au 26 


-0.02267(23) 
-0.01161(15) 
-0.01105(13) 


-0.02191(58) 
-0.01118(31) 
-0.01076(34) 


-0.02290(38) 
-0.01190(28) 
-0.01093(26) 


-0.02295(75) 
-0.01265(50) 
-0.01028(50) 


-0.0213(22) 
-0.0120(16) 
-0.00988(133) 


(u;2) 

<*% 
(d;2) 

av 26 
(2) 
av 26 


0.01943(28) 
0.01015(19) 
0.00926(15) 


0.01406(68) 
0.00797(33) 
0.00624(44) 


0.01122(48) 
0.00672(37) 
0.00445(32) 


0.00728(90) 
0.00611(66) 
0.00115(64) 


0.00206(265) 
0.00485(201) 
-0.00261(160) 


P = Pi 


«) 

(d) 
V 2b 

V2b 


0.4097(42) 

0.1906(23) 
0.2191(27) 


0.4059(97) 
0.1828(46) 
0.2231(61) 


0.4012(91) 
0.1813(46) 
0.2200(59) 


0.385(16) 
0.1697(72) 
0.215(12) 


0.390(34) 

0.162(15) 
0.228(25) 


(u;l) 
av 2b 

(d;l) 
av 26 

(1) 
au 26 


-0.02238(38) 
-0.01156(24) 
-0.01079(21) 


-0.0212(11) 
-0.01044(57) 
-0.01088(67) 


-0.02219(70) 
-0.01189(48) 
-0.01022(46) 


-0.0209(15) 
-0.01251(91) 
-0.00844(101) 


-0.0252(40) 
-0.0113(24) 
-0.0142(29) 


(«;2) 

aU 2fe 
(d;2) 

av 2b 
(2) 
aU 26 


0.01898(44) 
0.01003(28) 
0.00895(23) 


0.0142(13) 
0.00731(63) 
0.00706(79) 


0.01042(78) 
0.00675(56) 
0.00366(55) 


0.00566(163) 
0.00639(111) 
-0.00046(124) 


0.00748(445) 
0.00358(326) 
0.00485(456) 



Table 10: The bare results for v 2 from eq. (|39|) for (3 = 6.4, c sw = 1.526. 
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K 


0.1320 


0.1324 


0.1333 


0.1338 


0.1342 


P = Pl 


(u) 

Jd) 
v 3 

V3 


0.0981(55) 

0.0435(26) 
0.0546(34) 


0.0973(52) 

0.0408(25) 
0.0566(31) 


0.0961(72) 

0.0378(36) 
0.0582(46) 


0.114(15) 
0.0407(61) 
0.0710(98) 


0.0955(125) 
0.0342(62) 
0.0606(85) 


(«;mi) 
v 3 

v 3 

(mi) 
^3 


0.00051(215) 
0.00007(129) 
0.00045(150) 


0.00165(221) 
-0.00140(145) 
0.00056(155) 


0.00193(373) 
-0.00208(232) 
-0.00009(298) 


0.00694(638) 
-0.00158(379) 
0.0103(57) 


-0.00015(818) 
-0.00160(557) 
0.00113(687) 


(«;ma) 
u 3 

„,(d;m 2 ) 
v 3 

v 3 


-0.00490(717) 
-0.00342(459) 
-0.00134(487) 


-0.0102(73) 
-0.0108(54) 
0.00002(522) 


-0.0070(132) 

-0.0118(100) 
0.0022(112) 


0.0350(269) 
0.0053(167) 
0.0308(212) 


-0.0205(311) 
0.0146(227) 
0.0073(273) 




«) 
^4 


0.0272(28) 
0.0128(15) 
0.0143(19) 


0.0330(27) 
0.0128(16) 
0.0201(18) 


0.0342(39) 
0.0112(24) 
0.0232(27) 


0.0339(74) 
0.0141(42) 
0.0206(55) 


0.0331(64) 
0.00780(485) 
0.0253(47) 


(u:mi ) 

(d;mi) 
y 4 

(mi) 
V 4 


-0.00588(626) 
0.00355(394) 
-0.00905(399) 


-0.00625(648) 
-0.00113(450) 
-0.00512(392) 


-0.0048(104) 
0.00144(708) 
-0.00612(799) 


-0.0109(217) 
0.0086(129) 
-0.0202(144) 


0.0196(243) 
0.0152(153) 
0.0047(214) 


,,(«;"*a) 
v 4 
(d;m 2 ) 

vt 2) 


-0.00680(598) 
0.00180(298) 
-0.00803(382) 


-0.00273(503) 
-0.00188(359) 
-0.00169(327) 


-0.00053(841) 
-0.00142(591) 
0.00038(692) 


-0.0052(198) 
0.0066(100) 
-0.0085(150) 


0.0163(201) 
-0.0068(150) 
0.0243(180) 


v^ mz) /a 
vf mz) /a 


0.0178(75) 
0.00244(508) 
0.0156(48) 


0.00521(918) 
0.00523(630) 
-0.00127(587) 


0.0082(157) 
0.0093(108) 
-0.0041(107) 


0.0028(330) 
-0.0111(188) 
0.0150(255) 


0.0378(447) 
0.0434(496) 
-0.0062(351) 



Table 11: The bare results for ^3 and v$ from eq. (|39|) for (3 = 6.0, c sw = 1.769. 
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K 


0.1333 


0.1339 


0.1344 


0.1349 


P = Pl 


(u) 
v 3 

v 3 
V3 


0.1087(72) 

0.0493(36) 
0.0594(41) 


0.1075(64) 
0.0467(30) 
0.0608(40) 


0.110(11) 
0.0478(53) 
0.0628(71) 


0.115(12) 
0.0455(54) 
0.0697(88) 


v 3 

(d;mi) 
v 3 

(mi) 
v 3 


-0.00141(292) 
-0.00233(195) 
0.00062(175) 


0.00278(316) 
-0.00023(187) 
0.00294(225) 


0.00002(551) 
-0.00364(328) 
0.00364(415) 


-0.00474(849) 
-0.00716(453) 
0.00319(699) 


(u;m 2 ) 
v 3 

(d;m 2 ) 
v 3 

(m,2) 
v 3 


-0.0011(101) 
-0.00088(704) 
0.00191(505) 


-0.0038(120) 
-0.00597(810) 
0.00098(765) 


0.0201(228) 
0.0136(159) 
0.0135(143) 


0.0449(378) 
0.0262(224) 
0.0223(288) 


P = Pl 


(u) 

V A 
V 4 


0.0333(41) 
0.0127(21) 
0.0205(31) 


0.0336(55) 
0.0129(30) 
0.0206(37) 


0.0331(69) 

0.00757(388) 
0.0238(62) 


0.0442(85) 
0.00803(561) 
0.0344(89) 


(ti:mi ) 
(d;mi) 
(mi) 


0.0212(116) 
0.0102(69) 
0.0118(73) 


-0.0113(147) 
-0.0117(90) 
0.0083(906) 


0.0339(256) 
0.0176(136) 
0.0172(180) 


0.0375(332) 
0.0088(189) 
0.0270(267) 


V A 
(d;m 2 ) 


0.0074(101) 
0.00098(587) 
0.00756(606) 


-0.0007(130) 
-0.00384(687) 
0.00218(904) 


0.0157(213) 
0.0082(121) 
0.0100(150) 


0.0357(310) 
0.0092(168) 
0.0257(249) 


4" ;m3) /a 
vt m&) /a 
vt 3) /a 


0.0225(166) 
0.0090(110) 
0.0136(104) 


-0.0201(192) 
-0.00002(1285) 
-0.0134(129) 


0.0324(364) 
0.0198(235) 
0.0097(275) 


-0.0271(578) 
0.0291(390) 
-0.0520(497) 



Table 12: The bare results for ^3 and 1)4 from eq. (|39|) for (3 = 6.2, c sw = 1.614. 
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K 


0.1338 


0.1342 


0.1346 


0.1350 


0.1353 






P 


'=Pi 






(u) 
v 3 

v 3 
V3 


0.0996(62) 
0.0453(31) 
0.0543(37) 


0.0940(88) 
0.0413(43) 
0.0516(53) 


0.0963(88) 
0.0410(40) 
0.0554(59) 


0.0946(104) 
0.0381(54) 
0.0564(73) 


0.0898(179) 
0.0313(71) 
0.0594(148) 


\U'.Tn~\ ) 

v 3 

(d;mi) 
v 3 

Ami) 
v 3 


-0.00088(245) 
-0.00166(157) 
-0.00056(139) 


-0.00250(353) 
-0.00301(268) 
0.00018(258) 


-0.00078(418) 
-0.00270(278) 
0.00187(271) 


-0.00659(559) 
-0.00400(400) 
-0.00341(430) 


0.0180(117) 
-0.00394(813) 
0.0213(109) 


(u;m2) 
V 3 

v 3 

("12 ) 

v 3 


-0.0194(96) 
-0.0121(58) 
-0.00752(523) 


0.0135(151) 
-0.0016(110) 
0.0131(103) 


-0.0240(180) 
-0.0152(110) 
-0.0057(109) 


-0.0067(282) 
-0.0039(183) 
0.0014(191) 


0.0557(718) 
0.0212(400) 
0.0383(582) 






P 


'=Pi 






(«) 

v 4 
V4 


0.0355(49) 
0.0136(32) 
0.0216(27) 


0.0382(78) 
0.0136(37) 
0.0249(51) 


0.0440(75) 
0.0161(52) 
0.0280(46) 


0.0496(104) 
0.0169(72) 
0.0310(69) 


0.0403(227) 
0.0303(137) 
0.0123(166) 


lu:m\ ) 

v\ 
ld;mi) 

(mi) 

v\ 


-0.0067(130) 
-0.00681(859) 
-0.00076(782) 


0.0220(222) 
0.0162(132) 
-0.0006(131) 


-0.0004(217) 
-0.0117(151) 
0.0085(154) 


0.0194(296) 
0.0046(199) 
0.0107(230) 


0.0878(792) 
0.0106(445) 
0.0727(662) 


(u;m 2 ) 
v 4 

°A 


-0.0055(112) 
-0.00045(699) 
-0.00455(654) 


0.0189(190) 
0.0163(141) 
0.00174(936) 


0.0064(187) 
0.0050(117) 
-0.0001(127) 


0.0356(250) 
0.0264(164) 
0.0075(208) 


0.0102(674) 
0.0029(425) 
0.0038(508) 


vt m ' i] /a 
vt 3) /a 


0.0306(223) 
0.0182(124) 
0.0072(124) 


-0.0381(328) 
-0.0155(202) 
-0.0211(255) 


0.0570(402) 
0.0325(221) 
0.0100(259) 


0.142(57) 
0.0961(364) 
0.0488(449) 


-0.069(127) 
-0.100(84) 
0.028(108) 



Table 13: The bare results for ^3 and V4 from eq. (|3T))) for (3 = 6.4, c sw = 1.526. 
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